text:
import ProximalBundleMethod
import LinearAlgebra
import Random
using Printf

include("utils.jl")
include("csv_exporter.jl")
include("plot_utils.jl")

function least_squares_objective(x, A, b)
    err = A * x - b
    return LinearAlgebra.norm(err)
end

function least_squares_gradient(x, A, b)
    err = A * x - b
    if LinearAlgebra.norm(err) == 0.0
        return LinearAlgebra.zeros(length(err))
    end
    return A' * err / LinearAlgebra.norm(err)
end

function main()
    Random.seed!(926)
    n, m = 100, 60
    A = randn(n, m) / sqrt(n)
    x_opt = randn(m)
    b = A * x_opt
    x_init = zeros(m)

    objective = (x -> least_squares_objective(x, A, b))
    gradient = (x -> least_squares_gradient(x, A, b))
    step_size = (x, y, _) -> compute_ideal_step_size(x, y, 0.0, x_opt)

    params = create_bundle_method_parameters(
        iteration_limit = 150,
        verbose = true,
        printing_frequency = 50,
    )
    step_sizes = [1e-8 * 10^j for j = 0:8]

    println("About to solve a random least squares problem using ideal step size.")
    sol_ideal = ProximalBundleMethod.solve(objective, gradient, params, step_size, x_init)

    println(
        "\nAbout to solve a random least squares problem using adaptive parallel method.",
    )
    sol, iter_info_agents =
        ProximalBundleMethod.solve_adaptive(objective, gradient, params, step_sizes, x_init)

    csv_path_ideal = "results/regression/rusults_ideal_step_size_sharp_regression.csv"
    output_stepsize_ideal_plot = "results/regression/rusults_ideal_step_size_sharp_regression.pdf"
    output_objective_ideal_plot = "results/regression/rusults_ideal_objectives_sharp_regression.pdf"

    # Upper and lower bounds for the stepsize log plots
    stepsize_lower_bound = 0.8
    stepsize_upper_bound = 1.0e15

    # Upper and lower bounds for the objective function log plots
    obj_lower_bound = 1.5e-16
    obj_upper_bound = 2.0

    export_statistics(sol_ideal, csv_path_ideal)
    plot_step_sizes(
        csv_path_ideal,
        output_stepsize_ideal_plot;
        lower_bound = stepsize_lower_bound,
        upper_bound = stepsize_upper_bound,
    )
    plot_objective(
        csv_path_ideal,
        output_objective_ideal_plot;
        lower_bound = obj_lower_bound,
        upper_bound = obj_upper_bound,
    )

    csv_path = "results/regression/results_agents_sharp_regression.csv"
    output_path = "results/regression/results_agents_sharp_regression.pdf"
    export_losses_agents(
        iter_info_agents,
        [(@sprintf "%.0e" 1 / mu) for mu in step_sizes],
        csv_path,
    )
    plot_agent_losses(
        csv_path,
        output_path;
        lower_bound = obj_lower_bound,
        upper_bound = obj_upper_bound,
    )

    csv_path_ideal = "results/regression/rusults_parallel_bundle_sharp_regression.csv"
    output_stepsize_plot = "results/regression/rusults_parallel_stepsize_sharp_regression.pdf"
    output_objective_plot = "results/regression/rusults_parallel_objective_sharp_regression.pdf"

    export_statistics(sol, csv_path_ideal)
    plot_step_sizes(
        csv_path_ideal,
        output_stepsize_plot;
        lower_bound = stepsize_lower_bound,
        upper_bound = stepsize_upper_bound,
    )
    plot_objective(
        csv_path_ideal,
        output_objective_plot;
        lower_bound = obj_lower_bound,
        upper_bound = obj_upper_bound,
    )

end

main()


text:
using Printf

function main()
    xy_or_r = 2          # 1:x,y 2:r,t
    xnum = 200           #       2:rnum
    ynum = 100           #       2:tnum
    lenx = 1.0           #       2:inr
    leny = 2.0           #       2:outr
    st_ang  = pi/2       # 2:angle
    en_ang  = 0          # 2:angle

    outdir = "grid"
    make_dir(outdir)
    result = "result"
    make_dir(result)
    result = "post_result"
    make_dir(result)
    nodes, xnum_max, ynum_max = mk_gird(xnum,ynum,lenx,leny,xy_or_r,st_ang,en_ang,outdir)
    vecA(nodes,xnum_max,ynum_max,outdir)
end

function mk_gird(xnum,ynum,lenx,leny,xy_or_r,st_ang,en_ang,outdir)
    """
    nodes[i,j,k]
    i   : x,r方向の番号
    j   : y,theta方向の番号
    k=1 : 点のx座標
    k=2 : 点のy座標

    nodes[1,:]      : x方向境界
    nodes[:,1]      : y方向境界
    nodes[xnum+2,:] : x方向境界
    nodes[:,ynum+2] : y方向境界
    """
    xnum_max = xnum+1+2
    ynum_max = ynum+1+2
    nodes = zeros(xnum_max, ynum_max, 3)
    if xy_or_r == 1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                x = lenx/(xnum)*(i-2)
                y = leny/(ynum)*(j-2)
                nodes[i,j,1] = x
                nodes[i,j,2] = y
            end
        end
    elseif xy_or_r == 2
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                r = lenx + (leny-lenx)/ynum*(j-2)
                theta = st_ang - (st_ang-en_ang)/xnum*(i-2)
                x = r * cos(theta)
                y = r * sin(theta)
                nodes[i,j,1]=x
                nodes[i,j,2]=y
            end
        end
    end

    #=
    仮想セルの作成：現在あるセル境界線を延長して作成
    そのまま延長した際にセルが交差するのを防ぐため、仮想セルの大きさを100^-1倍にする.
    ベクトルで考えれば下記のような計算になる．たぶん
    また、[1,1]等の角にはダミー値が入っている.
    =#
    for j in 1:ynum_max
        for k in 1:2
            nodes[1,j,k]        =  1.01*nodes[2,j,k] - 0.01*nodes[3,j,k]
            nodes[xnum_max,j,k] =  1.01*nodes[xnum_max-1,j,k] - 0.01*nodes[xnum_max-2,j,k]
        end
    end
    for i in 1:xnum_max
        for k in 1:2
            nodes[i,1,k]        =  1.01*nodes[i,2,k] - 0.01*nodes[i,3,k]
            nodes[i,ynum_max,k] =  1.01*nodes[i,ynum_max-1,k] - 0.01*nodes[i,ynum_max-2,k]
        end
    end

    fff=outdir*"/nodes"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        for i in 1:xnum_max
            for j in 1:ynum_max
                x = @sprintf("%8.8e", nodes[i,j,1])
                y = @sprintf("%8.8e", nodes[i,j,2])
                write(f,string(i)*" "*string(j)*" "*x*" "*y*"\n")
            end
        end
    end
    println("write "*fff)

    # nodes_forvtk
    nodes_num = zeros(Int,xnum_max, ynum_max)
    fff=outdir*"/nodes_forvtk"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        k=1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                x = @sprintf("%8.8e", nodes[i,j,1])
                y = @sprintf("%8.8e", nodes[i,j,2])
                z = @sprintf("%8.8e", nodes[i,j,3])
                write(f,string(k)*" "*x*" "*y*" "*z*"\n")
                nodes_num[i,j] = k
                k = k+1
            end
        end
    end
    println("write "*fff)

    fff=outdir*"/nodesnum"
    open(fff,"w") do f
        write(f,"nodesnum: xnum_max, ynum_max\n")
        write(f,string(xnum_max)*" "*string(ynum_max)*"\n")
    end

    ### element ###
    fff=outdir*"/element_forvtk"
    open(fff,"w") do f
        write(f,"elements:cell_xnum, lup,rup,ldown,rdown \n")
        k=1
        for i in 2:xnum_max-2
            for j in 2:ynum_max-2
                d1 = @sprintf("%1.0f", nodes_num[i,j])
                d2 = @sprintf("%1.0f", nodes_num[i,j+1])
                d3 = @sprintf("%1.0f", nodes_num[i+1,j+1])
                d4 = @sprintf("%1.0f", nodes_num[i+1,j])
                write(f,string(k)*" "*d1*" "*d2*" "*d3*" "*d4*"\n")
                k = k+1
            end
        end
    end
    println("write "*fff)

    return  nodes,xnum_max,ynum_max
end

function vecA(nodes,xnum_max,ynum_max,outdir)

    vecAx=zeros(xnum_max, ynum_max-1, 2)
    for i in 1:xnum_max
        for j in 1:ynum_max-1
            # 2dim Ax=(y1-y3,x1-x3,0)
            x = nodes[i,j+1,2]-nodes[i,j,2]
            y = nodes[i,j,1]-nodes[i,j+1,1]
            vecAx[i,j,1] = x
            vecAx[i,j,2] = y
        end
    end

    vecAy=zeros(xnum_max-1, ynum_max, 2)
    for i in 1:xnum_max-1
        for j in 1:ynum_max
            # 2dim Ay=(y1-y2,x1-x2,0)
            x = nodes[i,j,2]-nodes[i+1,j,2]
            y = nodes[i+1,j,1]-nodes[i,j,1]
            vecAy[i,j,1] = x
            vecAy[i,j,2] = y
        end
    end

    fff=outdir*"/vecAx"
    open(fff,"w") do f
        write(f,"vecAx: xnum, ynum , x vec, y vec\n")
        for i in 1:xnum_max
            for j in 1:ynum_max-1
                x = @sprintf("%8.8f", vecAx[i,j,1])
                y = @sprintf("%8.8f", vecAx[i,j,2])
                write(f,string(i)*" "*string(j)*" "*x*" "*y*"\n")
            end
        end
    end
    println("write "*fff)

    fff=outdir*"/vecAy"
    open(fff,"w") do f
        write(f,"vecAy: xnum, ynum , x vec, y vec\n")
        for i in 1:xnum_max-1
            for j in 1:ynum_max
                x = @sprintf("%8.8f", vecAy[i,j,1])
                y = @sprintf("%8.8f", vecAy[i,j,2])
                write(f,string(i)*" "*string(j)*" "*x*" "*y*"\n")
            end
        end
    end
    println("write "*fff)
end

function make_dir(outdir)
    k = 0
    try rm(outdir,recursive=true)
    catch
        mkdir(outdir)
        k = 1
    end

    if k == 0
        mkdir(outdir)
    end
end

# ---------------------------------
main()

ext: jl


text:
function TimeSpanEstimatedTroposphericModel(arg0::EstimatedTroposphericModel)
    return TimeSpanEstimatedTroposphericModel((EstimatedTroposphericModel,), arg0)
end

function TimeSpanEstimatedTroposphericModel(arg0::EstimatedTroposphericModel, arg1::TimeScale)
    return TimeSpanEstimatedTroposphericModel((EstimatedTroposphericModel, TimeScale), arg0, arg1)
end

function add_tropospheric_model_valid_after(obj::TimeSpanEstimatedTroposphericModel, arg0::EstimatedTroposphericModel, arg1::AbsoluteDate)
    return jcall(obj, "addTroposphericModelValidAfter", void, (EstimatedTroposphericModel, AbsoluteDate), arg0, arg1)
end

function add_tropospheric_model_valid_before(obj::TimeSpanEstimatedTroposphericModel, arg0::EstimatedTroposphericModel, arg1::AbsoluteDate)
    return jcall(obj, "addTroposphericModelValidBefore", void, (EstimatedTroposphericModel, AbsoluteDate), arg0, arg1)
end

function compute_zenith_delay(obj::TimeSpanEstimatedTroposphericModel, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::FieldAbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2)
end

function compute_zenith_delay(obj::TimeSpanEstimatedTroposphericModel, arg0::jdouble, arg1::Vector{jdouble}, arg2::AbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{jdouble}, (jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2)
end

function extract_parameters(obj::TimeSpanEstimatedTroposphericModel, arg0::Vector{RealFieldElement}, arg1::FieldAbsoluteDate)
    return jcall(obj, "extractParameters", Vector{RealFieldElement}, (Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1)
end

function extract_parameters(obj::TimeSpanEstimatedTroposphericModel, arg0::Vector{jdouble}, arg1::AbsoluteDate)
    return jcall(obj, "extractParameters", Vector{jdouble}, (Vector{jdouble}, AbsoluteDate), arg0, arg1)
end

function get_parameters_drivers(obj::TimeSpanEstimatedTroposphericModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_transitions(obj::TimeSpanEstimatedTroposphericModel)
    return jcall(obj, "getTransitions", NavigableSet, ())
end

function get_tropospheric_model(obj::TimeSpanEstimatedTroposphericModel, arg0::AbsoluteDate)
    return jcall(obj, "getTroposphericModel", EstimatedTroposphericModel, (AbsoluteDate,), arg0)
end

function mapping_factors(obj::TimeSpanEstimatedTroposphericModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::TimeSpanEstimatedTroposphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::TimeSpanEstimatedTroposphericModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "pathDelay", RealFieldElement, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::TimeSpanEstimatedTroposphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "pathDelay", jdouble, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
include("../src/Opt.jl")
include("../src/Model.jl")
h = Model.hill(k=0.5,n=3)

@show Opt.find_min(h,[0.,1.])

@show Opt.find_max(hill,[0.,30.])

@show Opt.find_MinMax(hill, [0.,30.])

using Optim

xx(w) = w

function t(x;functor)

return functor(x)

Optim.optimize(h, 0.1, 0.2)

##
using Optim
using Test

struct y
    p
end

function (s::y)(x)
    return x*s.p
end

model = y(1.0)
model2(x;p=1.0) = x*p

@show  model(2.0)
@show  model2(2.0)

@test typeof(y) <: Any # Because Optim.jl accepts :Any object

Optim.optimize(model, 0.1, 0.2);
Optim.optimize(model2, 0.1, 0.2);


text:
using Test, CircStats

@testset "CircStats" begin
    av = [0:0.2:π;0.4π:0.05:0.6π;π:0.4:2π]
    am = [av av]
    aa = [am;;;am]
    wv = [0:0.1:0.5π;0.4π:0.05:0.6π;0:0.2:0.5π]
    wm = [wv wv]
    wa = [wm;;;wm]
    d = 0.1

    @test circ_r(av) ≈ 0.481125184946105
    @test circ_r(av;w=wv)  ≈ 0.634569454038054
    @test circ_r(av;w=wv,d)  ≈ 0.634833935115386

    circ_r(am)
    circ_r(am;w=wm)
    circ_r(am;w=wm,d)
    circ_r(am;w=wm,d,dims=2)

    circ_r(aa)
    circ_r(aa;w=wa)
    circ_r(aa;w=wa,d)
    circ_r(aa;w=wa,d,dims=3)

    μ, ul, ll = circ_mean(av)
    @test μ ≈ 1.568889360630885
    @test ul ≈ 2.037005335388447
    @test ll ≈ 1.100773385873324
    μ, ul, ll = circ_mean(av;w=wv)
    @test μ ≈ 1.691236041279729
    @test ul ≈ 2.018457532767024
    @test ll ≈ 1.364014549792434

    circ_mean(am)
    circ_mean(am;w=wm)

    circ_mean(aa)
    circ_mean(aa;w=wa)

    @test circ_median(av) ≈ 1.556637061435917

    s, s0 = circ_std(av)
    @test s ≈ 1.018699970603607
    @test s0 ≈ 1.209651009981690

    circ_std(am)
    circ_std(am;dims=2)
    circ_std(aa)
    circ_std(am;dims=3)

    S, s = circ_var(av)
    @test S ≈ 0.518874815053895
    @test s ≈ 1.037749630107790

    circ_var(am)
    circ_var(am;dims=2)
    circ_var(aa)
    circ_var(am;dims=3)

    α = [1,2,3];β=[2,3,4]
    circ_dist(α,π)
    circ_dist(α,β)

    circ_dist2(α)
    circ_dist2(α,β)

    mp, ρp, μp = circ_moment(av)
    @test mp ≈ 0.000917488892268 + 0.481124310135703im
    @test ρp ≈ 0.481125184946105
    @test μp ≈ 1.568889360630885
    mp, ρp, μp = circ_moment(av;cent=true)
    @test mp ≈ 0.481125184946105 - 0.000000000000000im
    @test ρp ≈ 0.481125184946105
    @test μp ≈ 3.741981750042832e-17

    mp, ρp, μp = circ_moment(am)
    mp, ρp, μp = circ_moment(am;cent=true)

    b, b0 = circ_skewness(av)
    @test b ≈ -0.005585405160200
    @test b0 ≈ -0.014943791328267
    b, b0 = circ_skewness(am)

    k, k0 = circ_kurtosis(av)
    @test k ≈ 0.315477455279792
    @test k0 ≈ 0.972747287142949
    k, k0 = circ_kurtosis(am)

    circ_stats(av)

    p,z = circ_rtest(av)
    @test p ≈ 1.247327496731614e-04
    @test z ≈ 8.564813412808679

    p,m = circ_otest(av)
    @test p ≈ 0.003445833222941
    @test m ≈ 7

    idx = [fill(1,18);fill(2,19)]
    p,F = circ_wwtest(av,idx)
    @test p ≈ 0.489205605248093
    @test F ≈ 0.488523443700675

    idp = [repeat([1,2,3,4],outer=9);1]
    idq = [repeat([1,2],outer=18);1]
    p,s = circ_hktest(av,idp,idq)
    @test p[1] ≈ 0.988563843496855
    @test p[2] ≈ 0.865090769456459
    @test isnan(p[3])
    @test s[1] ≈ 0.917000658130361
    @test s[2] ≈ 0.289841683535636
    @test isnan(s[3])

    p,f = circ_ktest(av,reverse(av))
    @test p ≈ 0.999999999999993
    @test f ≈ 1

    p = circ_symtest(av)
    # @test p ≈ 0.850454511827570 # diff in signrank in matlab?

    p,k,K=circ_kuipertest(av,reverse(av))
    @test p ≈ 1
    @test k ≈ 0

    p, u, UC = circ_raotest(av)
    @test p ≈ 0.5
    @test u ≈ 121.1858500639162
    @test UC ≈ 153.82

    p, v = circ_vtest(av)
    @test p ≈ 0.496851365629160
    @test v ≈ 0.033947089013911

    p = circ_medtest(av)
    @test p ≈ .0004719867429230360

    h, μ, ul, ll = circ_mtest(av)
    @test h
    @test μ ≈ 1.568889360630885
    @test ul ≈ 2.037005335388447
    @test ll ≈ 1.100773385873324

    ρ,p = circ_corrcc(av,reverse(av))
    @test ρ ≈ 0.204639214081096
    @test p ≈ 0.187332538620452

    ρ,p = circ_corrcl(av,reverse(av))
    @test ρ ≈ 0.591325482630456
    @test p ≈ 0.001551058328371

    @test circ_kappa(av) ≈ 1.095105628679724

    ϕ,c = circ_samplecdf(av,n=8)
    @test ϕ[4] ≈ 2.356194490192345
    @test c[4] ≈ 0.783783783783784

end


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, @get_package, @get_numbertype, @get_ndims, SUPPORTED_PACKAGES, PKG_CUDA, PKG_NONE, NUMBERTYPE_NONE, NDIMS_NONE
import ParallelStencil: @require, @symbols, longnameof
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. Reset of ParallelStencil" begin
            @testset "Reset if not initialized" begin
                @require !@is_initialized()
                @reset_parallel_stencil()
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
                @test @get_ndims() == $NDIMS_NONE
            end;
            @testset "Reset if initialized" begin
                @require !@is_initialized()
                @init_parallel_stencil($package, Float64, 3)
                @require @is_initialized() && @get_package() == $package
                @reset_parallel_stencil()
                @test length(@symbols($(@__MODULE__), Data)) == 1
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
                @test @get_ndims() == $NDIMS_NONE
            end;
        end;
    end;
)) end == nothing || true;


text:
"""
    bclogit(Y::Matrix{Int}, X::Matrix{Float64}, cs::Vector{Float64}; B::Int=2000, burn::Int=10000, printFreq::Int=0)

Baseline-category logit for Multinomial regression. Check out my notes at http://nbviewer.jupyter.org/github/luiarthur/GLM_AMS274/blob/master/notes/notes11.ipynb
"""
function bclogit(Y::Matrix{Int}, X::Matrix{Float64}, cs::Vector{Float64};
             B::Int=2000, burn::Int=10000, printFreq::Int=0)

  const (N,K) = size(X)

  assert( N == size(Y,1) )

  const J = size(Y,2)
  const stepSize_β = [cs[j] * sym(inv(X'X)) for j in 1:J-1]

  # can this computation be sped up?
  function ll(β::Matrix{Float64})
    EXB = exp(X * β)
    pi_denom = 1 + sum(EXB,2) # N by 1
    pi_mat = [ EXB ones(N) ] ./ pi_denom # N by J
    sum( Distributions.logpdf(Distributions.Multinomial(sum(Y[i,:]),pi_mat[i,:]), Y[i,:]) for i in 1:N )
  end

  function update(β::Matrix{Float64})
    newβ = copy(β)

    for j in 1:size(β,2)
      const cand = rand(Distributions.MvNormal(β[:,j], stepSize_β[j]))
      newβ[:,j] = cand

      if ll(newβ) - ll(β) < log(rand())
        newβ[:,j] = β[:,j]
      end
   end

   return newβ
  end

  const init = zeros(K,J-1)

  MCMC.gibbs(init, update, B, burn, printFreq=printFreq)

end


text:
using IterativeSolvers, KrylovKit, Arpack, LinearAlgebra

abstract type AbstractEigenSolver end
# abstract type for Matrix-Free eigensolvers
abstract type AbstractMFEigenSolver <: AbstractEigenSolver end
abstract type AbstractFloquetSolver <: AbstractEigenSolver end

# the following function returns the n-th eigenvectors computed by an eigen solver. This function is necessary given the different return types each eigensolver has
geteigenvector(eigsolve::ES, vecs, n::Union{Int, Array{Int64,1}}) where {ES <: AbstractEigenSolver} = vecs[:, n]
getsolver(eig::AbstractEigenSolver) = eig
####################################################################################################
# Solvers for default \ operator (backslash)
####################################################################################################
"""
The struct `Default` is used to  provide the backslash operator to our Package
"""
@with_kw struct DefaultEig{T} <: AbstractEigenSolver
	which::T = real		# how do we sort the computed eigenvalues
end

function (l::DefaultEig)(J, nev; kwargs...)
	# I put Array so we can call it on small sparse matrices
	F = eigen(Array(J))
	I = sortperm(F.values, by = l.which, rev = true)
	nev2 = min(nev, length(I))
	# we perform a conversion to Complex numbers here as the type can change from Float to Complex along the branch, this would cause a bug
	return Complex.(F.values[I[1:nev2]]), Complex.(F.vectors[:, I[1:nev2]]), true, 1
end

# case of sparse matrices or matrix free method
"""
$(TYPEDEF)
$(TYPEDFIELDS)

More information is available at [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl). You can pass the following parameters `tol=0.0, maxiter=300, ritzvec=true, v0=zeros((0,))`.

# Constructor

`EigArpack(sigma = nothing, which = :LR; kwargs...)`
"""
struct EigArpack{T, Tby, Tw} <: AbstractEigenSolver
	"Shift for Shift-Invert method with `(J - sigma⋅I)"
	sigma::T

	"Which eigen-element to extract :LR, :LM, ..."
	which::Symbol

	"Sorting function, default to real"
	by::Tby

	"Keyword arguments passed to EigArpack"
	kwargs::Tw
end

EigArpack(sigma = nothing, which = :LR; kwargs...) = EigArpack(sigma, which, real, kwargs)

function (l::EigArpack)(J, nev; kwargs...)
	if J isa AbstractMatrix
		λ, ϕ, ncv = Arpack.eigs(J; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	else
		if !(:v0 in keys(l.kwargs))
			error("The v0 argument must be provided in EigArpack for the matrix-free case")
		end
		N = length(l.kwargs[:v0])
		T = eltype(l.kwargs[:v0])
		Jmap = LinearMap{T}(J, N, N; ismutating = false)
		λ, ϕ, ncv, = Arpack.eigs(Jmap; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	end
	Ind = sortperm(λ, by = l.by, rev = true)
	ncv < nev && @warn "$ncv eigenvalues have converged using Arpack.eigs, you requested $nev"
	return λ[Ind], ϕ[:, Ind], true, 1
end
####################################################################################################
# Solvers for KrylovKit
####################################################################################################
"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
@with_kw struct EigKrylovKit{T, vectype} <: AbstractMFEigenSolver
	"Krylov Dimension"
	dim::Int64 = KrylovDefaults.krylovdim

	"Tolerance"
	tol::T = 1e-4

	"Number of restarts"
	restart::Int64 = 200

	"Maximum number of iterations"
	maxiter::Int64 = KrylovDefaults.maxiter

	"Verbosity ∈ {0,1,2}"
	verbose::Int = 0

	"Which eigenvalues are looked for :LR (largest real), :LM, ..."
	which::Symbol = :LR

	"If the linear map is symmetric, only meaningful if T<:Real"
	issymmetric::Bool = false

	"If the linear map is hermitian"
	ishermitian::Bool = false

	"Example of vector to usen for Krylov iterations"
	x₀::vectype = nothing
end

function (l::EigKrylovKit{T, vectype})(J, _nev; kwargs...) where {T, vectype}
	# note that there is no need to order the eigen-elements. KrylovKit does it
	# with the option `which`, by decreasing order.
	if J isa AbstractMatrix && isnothing(l.x₀)
		nev = min(_nev, size(J, 1))
		vals, vec, info = KrylovKit.eigsolve(J, nev, l.which;  verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, tol = l.tol, issymmetric = l.issymmetric, ishermitian = l.ishermitian)
	else
		nev = min(_nev, length(l.x₀))
		vals, vec, info = KrylovKit.eigsolve(J, l.x₀, nev, l.which;  verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, tol = l.tol, issymmetric = l.issymmetric, ishermitian = l.ishermitian)
	end
	# (length(vals) != _nev) && (@warn "EigKrylovKit returned $(length(vals)) eigenvalues instead of the $_nev requested")
	info.converged == 0 && (@warn "KrylovKit.eigsolve solver did not converge")
	return vals, vec, true, info.numops
end

geteigenvector(eigsolve::EigKrylovKit{T, vectype}, vecs, n::Union{Int, Array{Int64,1}}) where {T, vectype} = vecs[n]
####################################################################################################
# Solvers for ArnoldiMethod
####################################################################################################
# case of sparse matrices or matrix free method
"""
$(TYPEDEF)
$(TYPEDFIELDS)

More information is available at [ArnoldiMethod.jl](https://github.com/haampie/ArnoldiMethod.jl). For example, you can pass the parameters `tol, mindim, maxdim, restarts`.

# Constructor

`EigArnoldiMethod(;sigma = nothing, which = ArnoldiMethod.LR(), x₀ = nothing, kwargs...)`
"""
struct EigArnoldiMethod{T, Tby, Tw, Tkw, vectype} <: AbstractEigenSolver
	"Shift for Shift-Invert method"
	sigma::T

	"Which eigen-element to extract LR(), LM(), ..."
	which::Tw

	"how do we sort the computed eigenvalues, defaults to real"
	by::Tby

	"Key words arguments passed to EigArpack"
	kwargs::Tkw

	"Example of vector used for Krylov iterations"
	x₀::vectype
end

EigArnoldiMethod(;sigma = nothing, which = ArnoldiMethod.LR(), x₀ = nothing, kwargs...) = EigArnoldiMethod(sigma, which, real, kwargs, x₀)

function (l::EigArnoldiMethod)(J, nev; kwargs...)
	if J isa AbstractMatrix
		if isnothing(l.sigma)
			decomp, history = ArnoldiMethod.partialschur(J; nev = nev, which = l.which, l.kwargs...)
		else
			F = factorize(l.sigma * LinearAlgebra.I - J)
			Jmap = LinearMap{eltype(J)}((y, x) -> ldiv!(y, F, x), size(J, 1), ismutating=true)
			decomp, history = ArnoldiMethod.partialschur(Jmap; nev = nev, which = l.which, l.kwargs...)
		end
	else
		N = length(l.x₀)
		T = eltype(l.x₀)
		isnothing(l.sigma) == false && @warn "Shift-Invert strategy not implemented for maps"
		Jmap = LinearMap{T}(J, N, N; ismutating = false)
		decomp, history = ArnoldiMethod.partialschur(Jmap; nev = nev, which = l.which, l.kwargs...)
	end
	λ, ϕ = partialeigen(decomp)
	# shift and invert
	if isnothing(l.sigma) == false
		λ .= l.sigma .- 1 ./ λ
	end
	Ind = sortperm(λ, by = l.by, rev = true)
	ncv = length(λ)
	ncv < nev && @warn "$ncv eigenvalues have converged using ArnoldiMethod.partialschur, you requested $nev"
	return Complex.(λ[Ind]), Complex.(ϕ[:, Ind]), true, 1
end


text:
using Revise
using ADCME
using ADCMEKit
using NNFEM
using PyPlot
using ProgressMeter
using Printf
close("all")
include("load_domain_function.jl")

NT = 100
Δt = 20/NT

#################### solve static elastic response  ####################
domain = load_crack_domain()
H = domain.elements[1].mat[1].H
EBC_func = nothing
FBC_func = nothing
body_func = nothing
globaldata = GlobalData(missing, missing, missing, missing, domain.neqs, EBC_func, FBC_func, body_func)
Fext = compute_external_force(globaldata, domain)
d = LinearStaticSolver(globaldata, domain, domain.state, H, Fext)
sess = Session(); init(sess)
domain.state = run(sess, d)

figure()
subplot(211)
visualize_scalar_on_scoped_body(domain.state[1:domain.nnodes], zeros(size(domain.state)...), domain)
autoscale(enable=true, axis="both", tight=true)
xticks([])
xlabel("")
title(L"$u_x$")
subplot(212)
visualize_scalar_on_scoped_body(domain.state[domain.nnodes+1:end], zeros(size(domain.state)...), domain)
autoscale(enable=true, axis="both", tight=true)
title(L"$u_y$")
savefig("figures/dippingfault_disp_elasticity.png")

#################### solve viscoelastic relaxation  ####################
strain = getStrain(domain)
stress = getStress(domain)
state = domain.state

domain = load_crack_domain(option="mixed")
ts = GeneralizedAlphaSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)

d0 = state
v0 = zeros(2domain.nnodes)
a0  = zeros(2domain.nnodes)
# σ0 = zeros(getNGauss(domain),3)
# ϵ0 = zeros(getNGauss(domain),3)
σ0 = stress
ϵ0 = strain

μ = zeros(getNGauss(domain))
λ = zeros(getNGauss(domain))
η = zeros(getNGauss(domain))
k = 0
for i = 1:domain.neles
  e = domain.elements[i]
  for mat in e.mat
    global k += 1
    μ[k] = mat.μ
    λ[k] = mat.λ
    η[k] = mat.η
  end
end

d, v, a, σ, ϵ = ViscoelasticitySolver(
  globaldata, domain, d0, v0, a0, σ0, ϵ0, Δt, NT, μ, λ, η, Fext, ubd, abd
)

sess = Session(); init(sess)

d_, σ_ = run(sess, [d, σ])
domain.history["state"] = []
domain.history["stress"] = []
for i = 1:NT+1
  push!(domain.history["state"], d_[i,:])
  push!(domain.history["stress"], σ_[i,:,:])
end

#################### visulalize results ####################

## get surface slip
slip_idx = findall(domain.nodes[:,2] .≈ 0.0)
coords = domain.nodes[slip_idx,1]
ii = sortperm(coords)
coords = coords[ii]
x_ = d_[:, slip_idx[ii]]
y_ = d_[:, slip_idx[ii] .+ domain.nnodes]

x0 = d0[slip_idx[ii]]
y0 = d0[slip_idx[ii] .+ domain.nnodes]

figure(figsize=(6,3))
subplot(211)
plot(coords, (x_' .- x0)[:,1:20:end])
title(L"u_x - u_{x0}")
autoscale(enable=true, axis="x", tight=true)
xticks([])
subplot(212)
plot(coords, x_[1:20:end, :]')
title(L"u_x")
autoscale(enable=true, axis="x", tight=true)
tight_layout()
xlabel("x")
savefig("figures/dippingfault_ux_viscosity.png")

figure(figsize=(6,3))
subplot(211)
plot(coords, -(y_' .- y0)[:,1:20:end])
title(L"u_y - u_{y0}")
autoscale(enable=true, axis="x", tight=true)
xticks([])
subplot(212)
plot(coords, -y_[1:20:end, :]')
title(L"u_y")
autoscale(enable=true, axis="x", tight=true)
tight_layout()
xlabel("x")
savefig("figures/dippingfault_uy_viscosity.png")

figure(figsize=(5,3))
subplot(211)
pl1, = plot(coords, x0, "o-", markersize = 3)
autoscale(enable=true, axis="x", tight=true)
t = title("time = 0")
xticks([])
ylabel(L"u_x")
subplot(212)
pl2, = plot(coords, y0, "o-", markersize = 3)
autoscale(enable=true, axis="x", tight=true)
xlabel("x")
ylabel(L"u_y")
tight_layout()
function update(i)
  pl1.set_data(coords, x_[i, :])
  pl2.set_data(coords, -y_[i, :])
  t.set_text("time = $(round(i*Δt, digits=1))")
end
p = animate(update, collect(1:2:NT))
saveanim(p, "figures/dippingfault_viscosity.gif")

## get surface stress
idx = []
k = 0
pts = getGaussPoints(domain)
xcoord = Float64[]
for e in domain.elements
    global k
    if sum(map(x->x≈0.0, e.coords[:,2]))≠2
        k += length(e.weights)
        continue
    end
    l = argmin(pts[k+1:k+length(e.weights), 2])
    push!(idx, k + l)
    push!(xcoord, pts[k+l, 1])
    k += length(e.weights)
end
ii = sortperm(xcoord)
xcoord = xcoord[ii]
idx = idx[ii]

σ_surface = σ_[:, idx, :]
s = zeros(size(σ_surface, 1), size(σ_surface, 2))
for i = 1:size(σ_surface, 1)
    for j = 1:size(σ_surface, 2)
        s[i,j] = NNFEM.postprocess_stress(σ_surface[i,j,:], "vonMises")
    end
end

figure(figsize=(8,3))
plot(s[1:20:end, :]')
title("von Mises stress")
savefig("figures/dippingfault_vonMises_stress.png")

figure()
visualize_boundary(domain)
# visualize_mesh(domain)
visualize_scalar_on_scoped_body(η[1:9:end], zeros(domain.nnodes*2), domain)
gca().invert_yaxis()
# autoscale(enable=true, axis="both", tight=true)
scatter(domain.nodes[slip_idx,1], domain.nodes[slip_idx,2], color="m", marker=".", label="Tractrion Free")
legend()
title(L"Viscosity ($\eta$) model")
savefig("figures/dippingfault_mesh.png")

# figure()
# p = visualize_scalar_on_scoped_body(σ_[:, 1:domain.nnodes, 2], d_, domain, scale_factor=1.0)
# # p = visualize_von_mises_stress_on_scoped_body(d_, domain, scale_factor=10.0)
# # p = visualize_total_deformation_on_scoped_body(d_, domain, scale_factor=10.0)
# saveanim(p, "ad_solver.gif")
# close("all")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Iterators
using Distributions

export pred, pred_all
export solve_full

## supports probe_vec of DataFrame or Matrix{Integer}
function pred(r::Relation, probe_vec, F)
  if ! hasFeatures(r)
    return udot(r, probe_vec) + r.model.mean_value
  end
  return udot(r, probe_vec) + F * r.model.beta + r.model.mean_value
end

function pred(r::Relation)
  udot(r, r.data.df) + (hasFeatures(r) ? r.temp.linear_values : r.model.mean_value)
end

## computes predictions sum(u1 .* u2 .* ... .* uR, 1) for relation r
function udot(r::Relation, probe_vec::DataFrame)
  U = r.entities[1].model.sample[ :, probe_vec[:,1] ]
  for i in 2:length(r.entities)
    U .*= r.entities[i].model.sample[ :, probe_vec[:,i] ]
  end
  return vec(sum(U, 1))
end

## faster udot (about 2x) for matrices
function udot(r::Relation, probe_vec::Matrix)
  if length(size(r)) == 2
    ## special code for matrix
    num_latent = size(r.entities[1].model.sample, 1)
    result = zeros(size(probe_vec, 1))
    s1 = r.entities[1].model.sample
    s2 = r.entities[2].model.sample
    for i = 1:length(result)
      i1 = probe_vec[i,1]
      i2 = probe_vec[i,2]
      for k = 1:num_latent
        result[i] += s1[k,i1] * s2[k,i2]
      end
    end
    return result
  end
  U = r.entities[1].model.sample[:,probe_vec[:,1]]
  for i in 2:length(r.entities)
    U .*= r.entities[i].model.sample[:,probe_vec[:,i]]
  end
  return vec(sum(U, 1))
end

function udot2(uid::Vector, vid::Vector, U::Matrix, V::Matrix)
  size(U,1)   == size(V,1) || error("U and V must have the same number of rows.")
  length(uid) == length(vid) || error("uid and vid must have the same length.")

  num_latent = size(U,1)
  result     = zeros(eltype(U), length(uid))
  @inbounds for i = 1:length(uid)
    uu = uid[i]
    vv = vid[i]
    tmp = 0.0
    @inbounds @simd for k = 1:num_latent
      tmp += U[k,uu] * V[k,vv]
    end
    result[i] = tmp
  end
  return result
end

## computes predictions sum(u1 .* u2 .* ... .* uR, 2) for all points in relation r
function udot_all(r::Relation)
  ## matrix version:
  if length(r.entities) == 2
    return r.entities[1].model.sample' * r.entities[2].model.sample
  end

  ## TODO: make tensor version faster:
  U = zeros( Int64[en.count for en in r.entities]... )
  num_latent = length(r.entities[1].model.mu)
  for p in product( map(en -> 1:en.count, r.entities)... )
    x = ones(num_latent)
    for i in 1:length(p)
      x .*= r.entities[i].model.sample[:, p[i]]
    end
    U[p...] = sum(x)
  end
  return U
end

## predict all
function pred_all(r::Relation)
  if hasFeatures(r)
    error("Prediction of all elements is not possible when Relation has features.")
  end
  udot_all(r) + r.model.mean_value
end

function makeClamped(x, clamp::Vector{Float64})
  isempty(clamp) && return x

  x2 = copy(x)
  x2[x2 .< clamp[1]] = clamp[1]
  x2[x2 .> clamp[2]] = clamp[2]
  return x2
end

function clamp!(x, clamp::Vector{Float64})
  if ! isempty(clamp)
    x[x .< clamp[1]] = clamp[1]
    x[x .> clamp[2]] = clamp[2]
  end
  return x
end

function ConditionalNormalWishart(U::Matrix{Float64}, mu::Vector{Float64}, beta_0::Real, Tinv::Matrix{Float64}, nu::Real)
  N  = size(U, 2)
  NU = sum(U, 2)
  NS = U * U'

  nu_N   = nu + N
  beta_N = beta_0 + N
  mu_N   = (beta_0*mu + NU) / (beta_0 + N)
  T_N    = inv( Tinv + NS + beta_0 * mu * mu' - beta_N * mu_N * mu_N')

  NormalWishart(vec(mu_N), beta_N, T_N, nu_N)
end

function sample_alpha(alpha_lambda0::Float64, alpha_nu0::Float64, err::Vector{Float64})
  Λ  = alpha_lambda0 * eye(1)
  n  = length(err)
  SW = inv(inv(Λ) + err' * err)
  return rand(Wishart(alpha_nu0 + n, SW))[1]
end

function sample_lambda_beta(β::Matrix{Float64}, Lambda_u::Matrix{Float64}, ν::Float64, μ::Float64)
  νx = ν + size(β, 1) * size(β, 2)
  μx = μ * νx / (ν + μ * trace( (β'*β) * Lambda_u) )
  b  = νx / 2
  c  = 2*μx / νx
  return rand(Gamma(b, c))
end

function grab_col{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, col::Integer)
  r = A.colptr[col]:A.colptr[col+1]-1
  A.rowval[r], A.nzval[r]
end

function sample_latent_all2!(rel::Relation, dataRefs::Vector, procs::Vector{Int}, mode::Int,mu_u, Lambda_u::Matrix{Float64})
  Nprocs = length(procs)
  length(dataRefs) == Nprocs || error("Number of procs ($(Nprocs)) must equal the number of dataRefs ($(length(dataRefs))).")
  ## 1. create split ranges 1:length(procs):idmax
  sample_u = rel.entities[mode].model.sample
  ranges   = StepRange[ i:Nprocs:size(sample_u,2) for i in 1:Nprocs ]
  sample_m = Matrix{Float64}[en.model.sample for en in rel.entities]
  sample_m[mode] = zeros(0,0)

  mu_vector = typeof(mu_u) <: Vector

  ## 2. call sample_latent_range_ref at each worker with its dataRef
  @sync begin
    for i in 1:Nprocs
      @async begin
        if mu_vector
          sample_u[:, ranges[i]] = remotecall_fetch(procs[i], sample_latent_range_ref, ranges[i], dataRefs[i], mode, rel.model.mean_value, sample_m, rel.model.alpha, mu_u, Lambda_u)
        else
          sample_u[:, ranges[i]] = remotecall_fetch(procs[i], sample_latent_range_ref, ranges[i], dataRefs[i], mode, rel.model.mean_value, sample_m, rel.model.alpha, mu_u[:,ranges[i]], Lambda_u)
        end
      end
    end
  end
end

function sample_latent_range_ref(urange, Au_ref::RemoteRef, mode::Int, mean_rating, sample_mt, alpha, mu_u, Lambda_u)
  sample_m = length(sample_mt)==2 ? sample_mt[3 - mode] : sample_mt
  return sample_latent_range(urange, fetch(Au_ref), mode, mean_rating, sample_m, alpha, mu_u, Lambda_u)
end

## Sampling U, V for 2-way relation. Used by parallel code
## mu_u is Matrix of size num_latent x length(urange)
function sample_latent_range(urange, Au::FastIDF, mode::Int, mean_rating, sample_mt, alpha, mu_u::Matrix{Float64}, Lambda_u)
  U = zeros(size(mu_u, 1), length(urange))
  for i in 1:length(urange)
    u = urange[i]
    U[:,i] = sample_user_basic(u, Au, mode, mean_rating, sample_mt, alpha, mu_u[:,i], Lambda_u)
  end
  return U
end

## Sampling U, V for 2-way relation. Used by parallel code
function sample_latent_range(urange, Au::FastIDF, mode::Int, mean_rating, sample_mt, alpha, mu_u::Vector{Float64}, Lambda_u)
  U = zeros(length(mu_u), length(urange))
  for i in 1:length(urange)
    u = urange[i]
    U[:,i] = sample_user_basic(u, Au, mode, mean_rating, sample_mt, alpha, mu_u, Lambda_u)
  end
  return U
end

function sample_user_basic(uu::Integer, Au::FastIDF, mode::Int, mean_rating, sample_mt::Matrix{Float64}, alpha::Float64, mu_u::Vector{Float64}, Lambda_u::Matrix{Float64})
  id, v = getData(Au, mode, uu)
  ff = id[:, mode == 1 ? 2 : 1]

  rr = v - mean_rating
  MM = sample_mt[:, ff]

  covar = inv(Lambda_u + alpha * (MM*MM'))
  mu    = covar * (alpha * MM * rr + Lambda_u * mu_u)

  # Sample from normal distribution
  chol(covar)' * randn(length(mu_u)) + mu
end

## for Tensors
function sample_user_basic(uu::Integer, Au::FastIDF, mode::Int, mean_rating, sample_mt::Vector{Matrix{Float64}}, alpha::Float64, mu_u::Vector{Float64}, Lambda_u::Matrix{Float64})
  id, v = getData(Au, mode, uu)

  rr = v - mean_rating

  #ff = id[:, mode == 1 ? 2 : 1]
  #MM = sample_mt[:, ff]
  omodes = filter(i-> i!=mode, 1:length(sample_mt))
  j  = omodes[1]
  MM = sample_mt[j][ :, id[:,j] ]
  for j = omodes[2:end]
    MM .*= sample_mt[j][ :, id[:,j] ]
  end

  covar = inv(Lambda_u + alpha * (MM*MM'))
  mu    = covar * (alpha * MM * rr + Lambda_u * mu_u)

  # Sample from normal distribution
  chol(covar)' * randn(length(mu_u)) + mu
end

type Block
  ux::Vector{Int}   ## ids of the latent variables
  vx::Vector{Int}   ## ids of the other side
  Yma::Matrix{Float64} ## Y values w/o mean, size: length(v_idx) x length(u_idx)
end

function sample_users_blocked(block::Block, sample_mt::Matrix{Float64}, alpha::Float64, mu_u::Vector{Float64}, Lambda_u::Matrix{Float64})
  MM    = sample_mt[:, block.vx]
  covar = inv(Lambda_u + alpha * (MM*MM'))
  mu    = covar * (alpha * MM * block.Yma .+ Lambda_u * mu_u)

  # Sample from normal distribution
  chol(covar)' * randn(length(mu_u), size(mu, 2)) + mu
end

function sample_user2_all!(s::Entity, modes::Vector{Int64}, modes_other::Vector{Vector{Int64}})
  msample = s.model.sample
  mu      = s.model.mu
  for mm = 1:s.count
    msample[:, mm] = sample_user2(s, mm, mu, modes, modes_other)
  end
end

function sample_user2_all!(s::Entity, mu_matrix::Matrix{Float64}, modes::Vector{Int64}, modes_other::Vector{Vector{Int64}})
  msample = s.model.sample
  for mm = 1:s.count
    msample[:, mm] = sample_user2(s, mm, mu_matrix[:,mm], modes, modes_other)
  end
end

function sample_user2(s::Entity, i::Int, mu_si::Vector{Float64}, modes::Vector{Int64}, modes_other::Vector{Vector{Int64}})
  Lambda_si = copy(s.model.Lambda)
  mux       = s.model.Lambda * mu_si

  for r = 1:length(s.relations)
    rel = s.relations[r]
    df  = getData(rel.data, modes[r], i)
    rr  = convert(Array, df[:,end]) - (hasFeatures(rel) ? rel.temp.linear_values[getI(rel.data, modes[r], i)] : rel.model.mean_value)
    modes_o1 = modes_other[r][1]
    modes_o2 = modes_other[r][2:end]

    MM = rel.entities[modes_o1].model.sample[ :, df[:,modes_o1] ]
    for j = modes_o2
      MM .*= rel.entities[j].model.sample[ :, df[:,j] ]
    end
    Lambda_si += rel.model.alpha * MM * MM'
    mux       += rel.model.alpha * MM * rr
  end
  covar = inv(Lambda_si)
  mu    = covar * mux

  # Sample from normal distribution
  chol(covar)' * randn(length(mu)) + mu
end

function sample_beta(entity, sample_u_c, Lambda_u, lambda_beta, use_ff::Bool, tol=NaN )
  D, N = size(sample_u_c)
  numF = size(entity.F, 2)
  if isnan(tol)  ## default tolerance
    tol = eps() * numF
  end

  mv = MultivariateNormal(zeros(D), inv(Lambda_u) )
  ## TODO: using Ft_mul_Bt will be faster
  Ft_y = Ft_mul_B(entity, sample_u_c' + rand(mv, N)') + sqrt(lambda_beta) * rand(mv, numF)'
  #Ft_y = At_mul_B(entity.F, sample_u_c + rand(mv, N)') + sqrt(lambda_beta) * rand(mv, numF)'

  if use_ff
    beta = solve_full(entity.FF, Ft_y, lambda_beta)
  elseif ! isempty(entity.Frefs)
    beta = solve_cg2(entity.Frefs, Ft_y, lambda_beta, tol=tol)
  else
    error("No CG solver if entity has no Frefs or FF.")
    #beta = solve_cg(entity.F, Ft_y, lambda_beta)
  end
  return beta, Ft_y
end

function solve_full(FF, rhs, lambda_beta)
  FFreg = copy(FF)
  for i = 1:size(FFreg, 2)
    FFreg[i,i] += lambda_beta
  end
  return FFreg \ rhs
end

function sample_beta_rel(r::Relation)
  N, F = size(r.F)
  alpha  = r.model.alpha
  lambda = r.model.lambda_beta

  res  = getValues(r.data) - udot(r, r.data.df) - r.model.mean_value
  aFt_y = alpha * (r.F' * (res + alpha^(-0.5) * randn(N))) + sqrt(lambda) * randn(F)

  if isdefined(r.temp, :FF)
    ## using FF to compute beta_r
    K = alpha * r.temp.FF + lambda * speye(F)
    return K \ aFt_y
  else
    error("conjugate gradient unimplemented for sampling relation beta")
  end
end

function update_latent_prior!(en::Entity, full_lambda_u::Bool)
  mj = en.model

  local U::Matrix{Float64}
  nu   = mj.nu0
  Tinv = mj.WI

  if hasFeatures(en)
    mj.uhat = F_mul_beta(en)'
    U = mj.sample - mj.uhat
    if full_lambda_u
      nu   += size(mj.beta, 1)
      Tinv += mj.beta' * mj.beta * en.lambda_beta
    end
  else
    U = mj.sample
  end

  mj.mu, mj.Lambda = rand( ConditionalNormalWishart(U, mj.mu0, mj.b0, Tinv, nu) )
  return nothing
end

function update_beta!(en::Entity, tol::Float64)
  if ! hasFeatures( en )
    return nothing
  end
  mj = en.model
  mj.beta, rhs = sample_beta(en, mj.sample .- mj.mu, mj.Lambda, en.lambda_beta, en.use_FF, tol)
  if en.lambda_beta_sample
    en.lambda_beta = sample_lambda_beta(mj.beta, mj.Lambda, en.nu, en.mu)
  end
end


text:
"""
    IrregularGrid

## Fields
- **gridvals**: The values of the grid.
- **extrapolation_bc**: The extrapolation condition.
"""
struct IrregularGrid{BC <: BoundaryCondition} <: InterpolationGrid
    gridvals::Vector{Float64}
    extrapolation_bc::Union{BC, Float64}
end

Base.minimum(g::IrregularGrid) = minimum(g.gridvals)
Base.maximum(g::IrregularGrid) = maximum(g.gridvals)
Base.min(g::IrregularGrid) = minimum(g.gridvals)
Base.max(g::IrregularGrid) = maximum(g.gridvals)

Base.range(g::IrregularGrid) =
    throw(MethodError("step is not defined for irregular grids"))

Base.step(g::IrregularGrid) =
    throw(MethodError("step is not defined for irregular grids"))


text:
function format_check(io::IOBuffer, x::PTree, s::State)
    if length(s.doc.format_skips) == 0
        print_notcode(io, x, s)
        return
    end

    skip = s.doc.format_skips[1]
    line_range = x.startline:x.endline

    if s.on && skip[1] in line_range && skip[2] in line_range
        # weird corner case where off and on toggle
        # are in the same comment block
        x.endline = skip[1]
        print_notcode(io, x, s)
        write(io, skip[3])
        x.startline = skip[2]
        x.endline = skip[2]
        print_notcode(io, x, s)
    elseif s.on && skip[1] in line_range
        x.endline = skip[1]
        print_notcode(io, x, s)
        write(io, skip[3])
        s.on = false
    elseif !s.on && skip[2] in line_range
        deleteat!(s.doc.format_skips, 1)
        s.on = true
        # change the startline, otherwise lines
        # prior to in the NOTCODE node prior to
        # "format: on" will be reprinted
        x.startline = skip[2]
        print_notcode(io, x, s)
        # previous NEWLINE node won't be printed
    else
        print_notcode(io, x, s)
    end
end

function print_leaf(io::IOBuffer, x::PTree, s::State)
    if x.typ === NOTCODE
        format_check(io, x, s)
    elseif x.typ === INLINECOMMENT
        print_inlinecomment(io, x, s)
    else
        s.on && write(io, x.val)
    end
    s.line_offset += length(x)
end

function print_tree(io::IOBuffer, x::PTree, s::State)
    if is_leaf(x)
        print_leaf(io, x, s)
        return
    end
    print_tree(io, x.nodes, s, x.indent)
end

function print_tree(io::IOBuffer, nodes::Vector{PTree}, s::State, indent::Int)
    ws = repeat(" ", max(indent, 0))
    for (i, n) in enumerate(nodes)
        if is_leaf(n)
            print_leaf(io, n, s)
        elseif n.typ === CSTParser.StringH
            print_stringh(io, n, s)
        else
            print_tree(io, n, s)
        end

        if n.typ === NEWLINE && s.on && i < length(nodes)
            if is_closer(nodes[i+1]) ||
               nodes[i+1].typ === CSTParser.Block || nodes[i+1].typ === CSTParser.Begin
                write(io, repeat(" ", max(nodes[i+1].indent, 0)))
                s.line_offset = nodes[i+1].indent
            elseif !skip_indent(nodes[i+1])
                write(io, ws)
                s.line_offset = indent
            end
        end
    end
end

function print_stringh(io::IOBuffer, x::PTree, s::State)
    # The indent of StringH is set to the the offset
    # of when the quote is first encountered in the source file.

    # This difference notes if there is a change due to nesting.
    diff = x.indent + 1 - s.line_offset
    # @info "" length(x) s.line_offset

    # The new indent for the string is index of when a character in
    # the multiline string is FIRST encountered in the source file - the above difference
    # +1 since the character is 1 space after the indent
    x.indent = max(x.nodes[1].indent + 1 - diff, 0)
    print_tree(io, x, s)
end

function print_notcode(io::IOBuffer, x::PTree, s::State)
    s.on || return
    for l = x.startline:x.endline
        ws, v = get(s.doc.comments, l, (0, "\n"))
        v == "" && continue
        if l == x.endline && v[end] == '\n'
            v = v[1:prevind(v, end)]
        end
        ws > 0 && write(io, repeat(" ", ws))
        write(io, v)
        if l != x.endline && v[end] != '\n'
            write(io, "\n")
        end
    end
end

function print_inlinecomment(io::IOBuffer, x::PTree, s::State)
    s.on || return
    ws, v = get(s.doc.comments, x.startline, (0, ""))
    isempty(v) && return
    v = v[end] == '\n' ? v[nextind(v, 1):prevind(v, end)] : v
    ws > 0 && write(io, repeat(" ", ws))
    write(io, v)
end


text:
# Tests for doctest = :fix
#
# DOCUMENTER_TEST_DEBUG= JULIA_DEBUG=all julia test/doctests/fix/tests.jl
#
isdefined(@__MODULE__, :TestUtilities) || (include("../../TestUtilities.jl"); using .TestUtilities)
module DocTestFixTest
using Documenter, Test
using ..TestUtilities: @quietly

# Type to reliably show() objects across Julia versions:
@eval Main begin
    struct ShowWrap
        s :: String
    end
    Base.show(io::IO, x::ShowWrap) = write(io, x.s)
    const DocTestFixArray_1234 = Main.ShowWrap("4×1×1 Array{Int64,3}:\n[:, :, 1] =\n 1\n 2\n 3\n 4")
    const DocTestFixArray_2468 = Main.ShowWrap("4×1×1 Array{Int64,3}:\n[:, :, 1] =\n 2\n 4\n 6\n 8")
end

# The version check is necessary due to a behaviour change in https://github.com/JuliaLang/julia/pull/32851
mktempdir_nocleanup(dir) = VERSION >= v"1.3.0-alpha.112" ? mktempdir(dir, cleanup = false) : mktempdir(dir)

function normalize_line_endings(filename)
    s = read(filename, String)
    return replace(s, "\r\n" => "\n")
end

function test_doctest_fix(dir)
    srcdir = mktempdir_nocleanup(dir)
    builddir = mktempdir_nocleanup(dir)
    @debug "Testing doctest = :fix" srcdir builddir

    # Pkg.add changes permission of files to read-only,
    # so instead of copying them we read + write.
    src_jl = joinpath(srcdir, "src.jl")
    index_md = joinpath(srcdir, "index.md")
    write(index_md, normalize_line_endings(joinpath(@__DIR__, "broken.md")))
    write(src_jl, normalize_line_endings(joinpath(@__DIR__, "broken.jl")))

    # fix up
    include(joinpath(srcdir, "src.jl")); @eval import .Foo
    @debug "Running doctest/fix doctests with doctest=:fix"
    @quietly makedocs(sitename="-", modules = [Foo], source = srcdir, build = builddir, doctest = :fix)

    # test that strict = true works
    include(joinpath(srcdir, "src.jl")); @eval import .Foo
    @debug "Running doctest/fix doctests with doctest=true"
    @quietly makedocs(sitename="-", modules = [Foo], source = srcdir, build = builddir, strict = true)

    # also test that we obtain the expected output
    @test normalize_line_endings(index_md) == normalize_line_endings(joinpath(@__DIR__, "fixed.md"))
    @test normalize_line_endings(src_jl) == normalize_line_endings(joinpath(@__DIR__, "fixed.jl"))
end

@testset "doctest fixing" begin
    if haskey(ENV, "DOCUMENTER_TEST_DEBUG")
        # in this mode the directories remain
        test_doctest_fix(mktempdir_nocleanup(@__DIR__))
    else
        mktempdir(test_doctest_fix, @__DIR__)
    end
end

end # module


text:
function addcol!(b, A, idx)
    rows = rowvals(A)
    vals = nonzeros(A)
    for i in nzrange(A, idx)
        b[rows[i]] += vals[i]
    end
    return nothing
end

function sumcol(A, idx)
    vals = nonzeros(A)
    acc = 0
    for i in nzrange(A, idx)
        acc += vals[i]
    end
    return acc
end

function sumrows(A)
    n = size(A, 1)
    acc = zeros(eltype(A), n)
    for j in axes(A, 2)
        addcol!(acc, A, j)
    end
    return acc
end


text:
import Images
using OCReract
using Test

pkg_path = abspath(joinpath(dirname(pathof(OCReract)), ".."))
path_to_test_img = "$(pkg_path)/test/files/testocr.png"

function test_run_tesseract()
    # TODO: read this expected text from a text file in test/files
    expected_text = "This is a lot of 12 point text to test the\nocr code and see if it works on all types\nof file format.\n\nThe quick brown dog jumped over the\nlazy fox. The quick brown dog jumped\nover the lazy fox. The quick brown dog\njumped over the lazy fox. The quick\nbrown dog jumped over the lazy fox."

    res = run_tesseract(path_to_test_img, "/tmp/res")

    @test res == true

    result_txt = ""

    open("/tmp/res.txt", "r") do f
        result_txt = read(f, String)
    end

    @test strip(expected_text) == strip(result_txt)

end

function test_run_and_get_output()
    # TODO: read this expected text from a text file in test/files
    expected_text = "This is a lot of 12 point text to test the\nocr code and see if it works on all types\nof file format.\n\nThe quick brown dog jumped over the\nlazy fox. The quick brown dog jumped\nover the lazy fox. The quick brown dog\njumped over the lazy fox. The quick\nbrown dog jumped over the lazy fox."

    result_txt = run_and_get_output(Images.load(path_to_test_img))

    @test strip(expected_text) == strip(result_txt)

end

function test_path_exist()
    res = run_tesseract("/tmp/not_existing_image.png", "/tmp/res")
    @test res == false
end

@testset "RunTesseract" begin
    # Normal run
    test_run_tesseract()
    # Run and get
    test_run_and_get_output()
    # Again to test warning about already existing output file
    test_run_and_get_output()
    # Test with non existent image
    test_path_exist()
end


text:
module Day07

using AdventOfCode2019
using Combinatorics

function day07(input::String = readInput(joinpath(@__DIR__, "..", "data", "day07.txt")))
    program = parse.(Int, split(input, ","))

    maxOut1 = maximum(run_amplifiers(copy(program), phases, 0) for phases in permutations([0, 1, 2, 3, 4]))
    maxOut2 = maximum(run_feedback_loop(copy(program), phases, 0) for phases in permutations([5, 6, 7, 8, 9]))
    return [maxOut1, maxOut2]
end

function run_amplifiers(program::Array{Int,1}, phases::Array{Int,1}, inputA::Int)
    channels = [Channel{Int}(10) for i = 1:6]
    for i = 1:5
        put!(channels[i], phases[i])
    end
    put!(channels[1], inputA)

    for i = 1:5
        @async AdventOfCode2019.IntCode.run_program!(copy(program), channels[i], channels[i+1], nothing)
    end
    result = take!(channels[6])
    for i = 1:6
        close(channels[i])
    end
    return result
end

function run_feedback_loop(program::Array{Int,1}, phases::Array{Int,1}, inputA::Int)
    channels = [Channel{Int}(10) for i = 1:5]
    status = [Channel{Bool}(1) for i = 1:5]
    for i = 1:5
        put!(channels[i], phases[i])
    end
    put!(channels[1], inputA)

    for i = 1:5
        @async AdventOfCode2019.IntCode.run_program!(copy(program), channels[i], channels[mod1(i+1, 5)], status[i])
    end
    for i = 1:5
        take!(status[i])
    end

    result = take!(channels[1])
    for i = 1:5
        close(channels[i])
        close(status[i])
    end
    return result
end

end # module


text:
# Compute the Number of Times a Pattern Appears in a Text

function count_pattern(text, pattern)
    k = length(pattern)
    count = 0
    for ix in 1:(length(text)-k+1)
        substring = text[ix:ix+k-1]
        if substring == pattern
            count += 1
        end
    end
    count
end

function main()
    text = readline()
    pattern = readline()
    result = count_pattern(text, pattern)
    println(result)
end

main()


text:
# using Turing
using Distributions
using DynamicPPL
using IRTracker

@model coinflip(y) = begin
    p ~ Beta(1, 1)
    N = length(y)
    for i = 1:N
        y[i] ~ Bernoulli(p)
    end
end

model = coinflip([1,1,0])
vi = VarInfo()
# model(vi)

# ⟨model⟩(⟨vi⟩, ()...) = nothing
#   @1: [Arg:§1:%1] model
#   @2: [Arg:§1:%2] vi
#   @3: [§1:%3] ⟨SampleFromPrior⟩(, ()...) = SampleFromPrior()
#     @1: [Arg:§1:%1] @3#1 = SampleFromPrior
#     @2: [Const:§1:%2]
#     @3: [§1:%3] new(@2) = SampleFromPrior()
#     @4: [§1:&1] return @3 = SampleFromPrior()
#   @4: [§1:%4] @1(@2, @3, ()...) = nothing
#     @1: [Arg:§1:%1] @4#1 = model
#     @2: [Arg:§1:%2] @4#2 = vi
#     @3: [Arg:§1:%3] @4#3 = SampleFromPrior()
#     @4: [§1:%4] ⟨DefaultContext⟩(, ()...) = DefaultContext()
#       @1: [Arg:§1:%1] @4#1 = DefaultContext
#       @2: [Const:§1:%2]
#       @3: [§1:%3] new(@2) = DefaultContext()
#       @4: [§1:&1] return @3 = DefaultContext()
#     @5: [§1:%5] @1(, (@2, @3, @4)...) = nothing
#       @1: [Arg:§1:%1] @5#1 = model
#       @2: [Arg:§1:%2] @5#2 = (vi, SampleFromPrior(), DefaultContext())
#       @3: [§1:%3] ⟨NamedTuple⟩(, ()...) = NamedTuple()
#       @4: [§1:%4] ⟨pairs⟩(@3, ()...) = Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}()
#       @5: [§1:%5] ⟨tuple⟩(@4, @1) = (Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), model)
#       @6: [§1:%6] ⟨Core._apply⟩(⟨DynamicPPL.#_#3⟩, @5, @2) = nothing
#       @7: [§1:&1] return @6 = nothing
#     @6: [§1:&1] return @5 = nothing
#   @5: [§1:&1] return @4 = nothing

trace = track(DepthLimitContext(2), model.f, vi, SampleFromPrior(), DefaultContext(), model);

# ⟨var"##inner_function#426#6"()⟩(⟨vi⟩, ⟨SampleFromPrior()⟩, ⟨DefaultContext()⟩, model) = nothing
#   @1: [Arg:§1:%1] var"##inner_function#426#6"()
#   @2: [Arg:§1:%2] vi
#   @3: [Arg:§1:%3] SampleFromPrior()
#   @4: [Arg:§1:%4] DefaultContext()
#   @5: [Arg:§1:%5] model
#   @6: [§1:%6] ⟨getproperty⟩(@5, ⟨:args⟩) = (y = [1, 1, 0],)
#   @7: [§1:%7] ⟨getproperty⟩(@6, ⟨:y⟩) = [1, 1, 0]
#   @8: [§1:%8] ⟨typeof⟩(@7) = Array{Int64,1}
#   @9: [Const:§1:%9] Type{#s23} where #s23<:Union{AbstractFloat, AbstractArray}
#   @10: [§1:%10] ⟨isa⟩(@7, @9) = false
#   @11: [§1:&1] goto §3 since @10 == false
#   @12: [Const:§3:%13] DynamicPPL.hasmissing
#   @13: [§3:%14] @12(@8) = false
#   @14: [§3:&1] goto §5 since @13 == false
#   @15: [§5:&1] goto §6 (@7)
#   @16: [Arg:§6:%18] @15#1 = [1, 1, 0]
#   @17: [Const:§6:%19] DynamicPPL.resetlogp!
#   @18: [§6:%20] @17(@2) = 0
#   @19: [§6:%21] ⟨Beta⟩(⟨1⟩, ⟨1⟩) = Beta{Float64}(α=1.0, β=1.0)
#   @20: [Const:§6:%22] (:msg,)
#   @21: [§6:%23] ⟨Core.apply_type⟩(⟨NamedTuple⟩, @20) = NamedTuple{(:msg,),T} where T<:Tuple
#   @22: [§6:%24] ⟨tuple⟩(⟨"Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340."⟩) = ("Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @23: [§6:%25] @21(@22) = (msg = "Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @24: [Const:§6:%26] DynamicPPL.assert_dist
#   @25: [§6:%27] ⟨Core.kwfunc⟩(@24) = DynamicPPL.var"#kw##assert_dist"()
#   @26: [Const:§6:%28] DynamicPPL.assert_dist
#   @27: [§6:%29] @25(@23, @26, @19) = true
#   @28: [§6:%30] ⟨Val⟩(⟨:p⟩) = Val{:p}()
#   @29: [Const:§6:%31] DynamicPPL.inparams
#   @30: [§6:%32] @29(@28, ⟨Val{(:y,)}()⟩) = false
#   @31: [§6:%33] ⟨!⟩(@30) = true
#   @32: [§6:&2] goto §7
#   @33: [§7:&1] goto §9 (@31)
#   @34: [Arg:§9:%38] @33#1 = true
#   @35: [§9:&2] goto §10
#   @36: [Const:§10:%39] VarName
#   @37: [§10:%40] ⟨Core.apply_type⟩(@36, ⟨:p⟩) = VarName{:p}
#   @38: [§10:%41] @37(⟨""⟩) = VarName{:p}("")
#   @39: [§10:%42] ⟨tuple⟩(@38, ⟨()⟩) = (VarName{:p}(""), ())
#   @40: [§10:&1] goto §18 (@39)
#   @41: [Arg:§18:%53] @40#1 = (VarName{:p}(""), ())
#   @42: [§18:%54] ⟨isa⟩(@41, ⟨Tuple⟩) = true
#   @43: [§18:&2] goto §19
#   @44: [§19:%55] ⟨Base.indexed_iterate⟩(@41, ⟨1⟩) = (VarName{:p}(""), 2)
#   @45: [§19:%56] ⟨getfield⟩(@44, ⟨1⟩) = VarName{:p}("")
#   @46: [§19:%57] ⟨getfield⟩(@44, ⟨2⟩) = 2
#   @47: [§19:%58] ⟨Base.indexed_iterate⟩(@41, ⟨2⟩, @46) = ((), 3)
#   @48: [§19:%59] ⟨getfield⟩(@47, ⟨1⟩) = ()
#   @49: [Const:§19:%60] DynamicPPL.tilde
#   @50: [§19:%61] @49(@4, @3, @19, @45, @48, @2) = (0.8047069487761294, -0.0)
#   @51: [§19:%62] ⟨getindex⟩(@50, ⟨1⟩) = 0.8047069487761294
#   @52: [Const:§19:%63] DynamicPPL.acclogp!
#   @53: [§19:%64] ⟨getindex⟩(@50, ⟨2⟩) = -0.0
#   @54: [§19:%65] @52(@2, @53) = 0.0
#   @55: [§19:&1] goto §21 (@51)
#   @56: [Arg:§21:%70] @55#1 = 0.8047069487761294
#   @57: [§21:%71] ⟨length⟩(@16) = 3
#   @58: [§21:%72] ⟨Colon()⟩(⟨1⟩, @57) = 1:3
#   @59: [§21:%73] ⟨iterate⟩(@58) = (1, 1)
#   @60: [§21:%74] ⟨===⟩(@59, ⟨nothing⟩) = false
#   @61: [§21:%75] ⟨not_int⟩(@60) = true
#   @62: [§21:&2] goto §22 (@59)
#   @63: [Arg:§22:%76] @62#1 = (1, 1)
#   @64: [§22:%77] ⟨getfield⟩(@63, ⟨1⟩) = 1
#   @65: [§22:%78] ⟨getfield⟩(@63, ⟨2⟩) = 1
#   @66: [§22:%79] ⟨Bernoulli⟩(@56) = Bernoulli{Float64}(p=0.8047069487761294)
#   @67: [Const:§22:%80] (:msg,)
#   @68: [§22:%81] ⟨Core.apply_type⟩(⟨NamedTuple⟩, @67) = NamedTuple{(:msg,),T} where T<:Tuple
#   @69: [§22:%82] ⟨tuple⟩(⟨"Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340."⟩) = ("Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @70: [§22:%83] @68(@69) = (msg = "Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @71: [Const:§22:%84] DynamicPPL.assert_dist
#   @72: [§22:%85] ⟨Core.kwfunc⟩(@71) = DynamicPPL.var"#kw##assert_dist"()
#   @73: [Const:§22:%86] DynamicPPL.assert_dist
#   @74: [§22:%87] @72(@70, @73, @66) = true
#   @75: [§22:%88] ⟨Val⟩(⟨:y⟩) = Val{:y}()
#   @76: [Const:§22:%89] DynamicPPL.inparams
#   @77: [§22:%90] @76(@75, ⟨Val{(:y,)}()⟩) = true
#   @78: [§22:%91] ⟨!⟩(@77) = false
#   @79: [§22:&1] goto §24 since @78 == false
#   @80: [Const:§24:%92] DynamicPPL.inparams
#   @81: [Const:§24:%93] DynamicPPL.getmissing
#   @82: [§24:%94] @81(@5) = Val{()}()
#   @83: [§24:%95] @80(@75, @82) = false
#   @84: [§24:&1] goto §25 (@83)
#   @85: [Arg:§25:%96] @84#1 = false
#   @86: [§25:&1] goto §30 since @85 == false
#   @87: [Const:§30:%116] DynamicPPL.inparams
#   @88: [§30:%117] @87(@75, ⟨Val{(:y,)}()⟩) = true
#   @89: [§30:&2] goto §31
#   @90: [§31:%118] ⟨getindex⟩(@16, @64) = 1
#   @91: [§31:%119] ⟨===⟩(@90, ⟨missing⟩) = false
#   @92: [§31:&1] goto §36 since @91 == false
#   @93: [§36:&1] goto §37 (@90)
#   @94: [Arg:§37:%139] @93#1 = 1
#   @95: [§37:&1] goto §39 (@94)
#   @96: [Arg:§39:%141] @95#1 = 1
#   @97: [§39:&1] goto §40 (@96)
#   @98: [Arg:§40:%142] @97#1 = 1
#   @99: [§40:%143] ⟨isa⟩(@98, ⟨Tuple⟩) = false
#   @100: [§40:&1] goto §42 since @99 == false
#   @101: [Const:§42:%156] DynamicPPL.acclogp!
#   @102: [Const:§42:%157] DynamicPPL.tilde
#   @103: [§42:%158] @102(@4, @3, @66, @98, @2) = -0.21727710662919572
#   @104: [§42:%159] @101(@2, @103) = -0.21727710662919572
#   @105: [§42:&1] goto §43
#   @106: [§43:%160] ⟨iterate⟩(@58, @65) = (2, 2)
#   @107: [§43:%161] ⟨===⟩(@106, ⟨nothing⟩) = false
#   @108: [§43:%162] ⟨not_int⟩(@107) = true
#   @109: [§43:&2] goto §22 (@106)
#   @110: [Arg:§22:%76] @109#1 = (2, 2)
#   @111: [§22:%77] ⟨getfield⟩(@110, ⟨1⟩) = 2
#   @112: [§22:%78] ⟨getfield⟩(@110, ⟨2⟩) = 2
#   @113: [§22:%79] ⟨Bernoulli⟩(@56) = Bernoulli{Float64}(p=0.8047069487761294)
#   @114: [Const:§22:%80] (:msg,)
#   @115: [§22:%81] ⟨Core.apply_type⟩(⟨NamedTuple⟩, @114) = NamedTuple{(:msg,),T} where T<:Tuple
#   @116: [§22:%82] ⟨tuple⟩(⟨"Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340."⟩) = ("Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @117: [§22:%83] @115(@116) = (msg = "Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @118: [Const:§22:%84] DynamicPPL.assert_dist
#   @119: [§22:%85] ⟨Core.kwfunc⟩(@118) = DynamicPPL.var"#kw##assert_dist"()
#   @120: [Const:§22:%86] DynamicPPL.assert_dist
#   @121: [§22:%87] @119(@117, @120, @113) = true
#   @122: [§22:%88] ⟨Val⟩(⟨:y⟩) = Val{:y}()
#   @123: [Const:§22:%89] DynamicPPL.inparams
#   @124: [§22:%90] @123(@122, ⟨Val{(:y,)}()⟩) = true
#   @125: [§22:%91] ⟨!⟩(@124) = false
#   @126: [§22:&1] goto §24 since @125 == false
#   @127: [Const:§24:%92] DynamicPPL.inparams
#   @128: [Const:§24:%93] DynamicPPL.getmissing
#   @129: [§24:%94] @128(@5) = Val{()}()
#   @130: [§24:%95] @127(@122, @129) = false
#   @131: [§24:&1] goto §25 (@130)
#   @132: [Arg:§25:%96] @131#1 = false
#   @133: [§25:&1] goto §30 since @132 == false
#   @134: [Const:§30:%116] DynamicPPL.inparams
#   @135: [§30:%117] @134(@122, ⟨Val{(:y,)}()⟩) = true
#   @136: [§30:&2] goto §31
#   @137: [§31:%118] ⟨getindex⟩(@16, @111) = 1
#   @138: [§31:%119] ⟨===⟩(@137, ⟨missing⟩) = false
#   @139: [§31:&1] goto §36 since @138 == false
#   @140: [§36:&1] goto §37 (@137)
#   @141: [Arg:§37:%139] @140#1 = 1
#   @142: [§37:&1] goto §39 (@141)
#   @143: [Arg:§39:%141] @142#1 = 1
#   @144: [§39:&1] goto §40 (@143)
#   @145: [Arg:§40:%142] @144#1 = 1
#   @146: [§40:%143] ⟨isa⟩(@145, ⟨Tuple⟩) = false
#   @147: [§40:&1] goto §42 since @146 == false
#   @148: [Const:§42:%156] DynamicPPL.acclogp!
#   @149: [Const:§42:%157] DynamicPPL.tilde
#   @150: [§42:%158] @149(@4, @3, @113, @145, @2) = -0.21727710662919572
#   @151: [§42:%159] @148(@2, @150) = -0.43455421325839144
#   @152: [§42:&1] goto §43
#   @153: [§43:%160] ⟨iterate⟩(@58, @112) = (3, 3)
#   @154: [§43:%161] ⟨===⟩(@153, ⟨nothing⟩) = false
#   @155: [§43:%162] ⟨not_int⟩(@154) = true
#   @156: [§43:&2] goto §22 (@153)
#   @157: [Arg:§22:%76] @156#1 = (3, 3)
#   @158: [§22:%77] ⟨getfield⟩(@157, ⟨1⟩) = 3
#   @159: [§22:%78] ⟨getfield⟩(@157, ⟨2⟩) = 3
#   @160: [§22:%79] ⟨Bernoulli⟩(@56) = Bernoulli{Float64}(p=0.8047069487761294)
#   @161: [Const:§22:%80] (:msg,)
#   @162: [§22:%81] ⟨Core.apply_type⟩(⟨NamedTuple⟩, @161) = NamedTuple{(:msg,),T} where T<:Tuple
#   @163: [§22:%82] ⟨tuple⟩(⟨"Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340."⟩) = ("Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @164: [§22:%83] @162(@163) = (msg = "Right-hand side of a ~ must be subtype of Distribution or a vector of Distributions on line 340.",)
#   @165: [Const:§22:%84] DynamicPPL.assert_dist
#   @166: [§22:%85] ⟨Core.kwfunc⟩(@165) = DynamicPPL.var"#kw##assert_dist"()
#   @167: [Const:§22:%86] DynamicPPL.assert_dist
#   @168: [§22:%87] @166(@164, @167, @160) = true
#   @169: [§22:%88] ⟨Val⟩(⟨:y⟩) = Val{:y}()
#   @170: [Const:§22:%89] DynamicPPL.inparams
#   @171: [§22:%90] @170(@169, ⟨Val{(:y,)}()⟩) = true
#   @172: [§22:%91] ⟨!⟩(@171) = false
#   @173: [§22:&1] goto §24 since @172 == false
#   @174: [Const:§24:%92] DynamicPPL.inparams
#   @175: [Const:§24:%93] DynamicPPL.getmissing
#   @176: [§24:%94] @175(@5) = Val{()}()
#   @177: [§24:%95] @174(@169, @176) = false
#   @178: [§24:&1] goto §25 (@177)
#   @179: [Arg:§25:%96] @178#1 = false
#   @180: [§25:&1] goto §30 since @179 == false
#   @181: [Const:§30:%116] DynamicPPL.inparams
#   @182: [§30:%117] @181(@169, ⟨Val{(:y,)}()⟩) = true
#   @183: [§30:&2] goto §31
#   @184: [§31:%118] ⟨getindex⟩(@16, @158) = 0
#   @185: [§31:%119] ⟨===⟩(@184, ⟨missing⟩) = false
#   @186: [§31:&1] goto §36 since @185 == false
#   @187: [§36:&1] goto §37 (@184)
#   @188: [Arg:§37:%139] @187#1 = 0
#   @189: [§37:&1] goto §39 (@188)
#   @190: [Arg:§39:%141] @189#1 = 0
#   @191: [§39:&1] goto §40 (@190)
#   @192: [Arg:§40:%142] @191#1 = 0
#   @193: [§40:%143] ⟨isa⟩(@192, ⟨Tuple⟩) = false
#   @194: [§40:&1] goto §42 since @193 == false
#   @195: [Const:§42:%156] DynamicPPL.acclogp!
#   @196: [Const:§42:%157] DynamicPPL.tilde
#   @197: [§42:%158] @196(@4, @3, @160, @192, @2) = -1.6332540217433904
#   @198: [§42:%159] @195(@2, @197) = -2.067808235001782
#   @199: [§42:&1] goto §43
#   @200: [§43:%160] ⟨iterate⟩(@58, @159) = nothing
#   @201: [§43:%161] ⟨===⟩(@200, ⟨nothing⟩) = true
#   @202: [§43:%162] ⟨not_int⟩(@201) = false
#   @203: [§43:&1] goto §44 since @202 == false
#   @204: [§44:&1] return ⟨nothing⟩


text:
"""
    Heatmap(xedges, yedges; left = true, closed = true)
    Heatmap(itr; left = true, closed = true)

Create a two dimensional histogram with the bin partition created by `xedges` and `yedges`.
When fitting a new observation, the first value will be associated with X, the second with Y.

- If `left`, the bins will be left-closed.
- If `closed`, the bins on the ends will be closed.  See [Hist](@ref).

# Example

    using Plots

    xy = zip(randn(10^6), randn(10^6))
    o = fit!(HeatMap(-5:.1:5, -5:.1:5), xy)
    plot(o)

    xy = zip(1 .+ randn(10^6) ./ 10, randn(10^6))
    o = HeatMap(xy)
    plot(o, marginals=false)
    plot(o)
"""
mutable struct HeatMap{EX, EY} <: OnlineStat{TwoThings}
    xedges::EX
    yedges::EY
    counts::Matrix{Int}
    out::Int
    left::Bool
    closed::Bool

    function HeatMap(x::T, y::S; left::Bool=true, closed::Bool=true) where {T<:AbstractVector,S<:AbstractVector}
        new{T,S}(x,y, zeros(Int, length(x)-1, length(y)-1), 0, left, closed)
    end
end
function HeatMap(itr, bins::Integer=200; kw...)
    xa, xb = extrema(first, itr)
    ya, yb = extrema(last, itr)
    fit!(HeatMap(range(xa, stop=xb, length=bins), range(ya, stop=yb, length=bins); kw...), itr)
end
nobs(o::HeatMap) = sum(o.counts) + o.out
value(o::HeatMap) = (x=o.xedges, y=o.yedges, z=o.counts)

function _fit!(o::HeatMap, xy)
    i = binindex(o.xedges, xy[1], o.left, o.closed)
    if 1 ≤ i < length(o.xedges)
        j = binindex(o.yedges, xy[2], o.left, o.closed)
        if 1 ≤ j < length(o.yedges)
            o.counts[i, j] += 1
        else
            o.out += 1
        end
    else
        o.out += 1
    end
end

function _merge!(a::HeatMap, b::HeatMap)
    if (a.xedges == b.xedges) && (a.yedges == b.yedges) && (a.left == b.left) && (a.closed == b.closed)
        a.out += b.out
        a.counts .+= b.counts
    else
        @warn "HeatMaps differ in edges or direction of bins being closed.  No merging occurred."
    end
end

@recipe function f(o::HeatMap; marginals=true)
    link --> :x
    z = Float64.(o.counts)
    z[z .== 0] .= NaN
    @series begin
        s = marginals ? 3 : 1
        subplot --> s
        seriestype --> :heatmap
        legend --> false
        o.xedges, o.yedges, z'
    end
    if marginals
        layout --> (2, 2)
        @series begin
            subplot --> 2
            seriestype --> :heatmap
            legendposition --> :left
            gride --> false
            axis --> false
            framestyle --> :none
            o.xedges, o.yedges, z'
        end
        @series begin
            subplot --> 1
            label --> ""
            linewidth --> 0
            line_alpha --> 0
            seriestype --> :bar
            o.xedges, vec(sum(o.counts, dims=1))
        end
        @series begin
            subplot --> 4
            label --> ""
            linewidth --> 0
            line_alpha --> 0
            orientation --> :h
            seriestype --> :bar
            o.yedges, vec(sum(o.counts, dims=2))
        end
    end
end


text:
# ---
# title: 72. Edit Distance
# id: problem72
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/edit-distance/description/>
# hidden: true
# ---
#
# Given two strings `word1` and `word2`, return _the minimum number of
# operations required to convert`word1` to `word2`_.
#
# You have the following three operations permitted on a word:
#
#   * Insert a character
#   * Delete a character
#   * Replace a character
#
#
#
# **Example 1:**
#
#
#
#     Input: word1 = "horse", word2 = "ros"
#     Output: 3
#     Explanation:
#     horse -> rorse (replace 'h' with 'r')
#     rorse -> rose (remove 'r')
#     rose -> ros (remove 'e')
#
#
# **Example 2:**
#
#
#
#     Input: word1 = "intention", word2 = "execution"
#     Output: 5
#     Explanation:
#     intention -> inention (remove 't')
#     inention -> enention (replace 'i' with 'e')
#     enention -> exention (replace 'n' with 'x')
#     exention -> exection (replace 'n' with 'c')
#     exection -> execution (insert 'u')
#
#
#
#
# **Constraints:**
#
#   * `0 <= word1.length, word2.length <= 500`
#   * `word1` and `word2` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm
import Base.show

function show(filename)
	V, EV = Plasm.svg2lar(filename)
	Plasm.view(V,EV)
	return V, EV
end

#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/new.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/curved.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/twopaths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/paths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/boundarytest2.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/tile.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/interior.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/holes.svg")

V,EV = show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/Lar.svg")


text:
module Algorithms

using Distributions: Categorical

using ..Scruff
using ..Utils
using ..RTUtils
using ..Models
using ..SFuncs
using ..Operators

include("algorithms/algorithm.jl")
include("algorithms/query.jl")
include("algorithms/sample_utils.jl")
include("algorithms/instant/instantalgorithm.jl")
include("algorithms/instant/ve.jl")
include("algorithms/instant/bp.jl")
include("algorithms/instant/threepassbp.jl")
include("algorithms/instant/loopybp.jl")
include("algorithms/instant/importance.jl")
include("algorithms/iterative/iterativealgorithm.jl")
include("algorithms/iterative/lazyinference.jl")
include("algorithms/iterative/lsfi.jl")
include("algorithms/iterative/iterativesampler.jl")
include("algorithms/filtering/filter.jl")
include("algorithms/filtering/windowutils.jl")
include("algorithms/filtering/windowcreator.jl")
include("algorithms/filtering/windowfilter.jl")
include("algorithms/filtering/bpfilter.jl")
include("algorithms/filtering/loopyfilter.jl")
include("algorithms/filtering/particlefilter.jl")
end


text:
module NumericalWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbsoluteJacobiansMapper = @jimport org.orekit.propagation.numerical.AbsoluteJacobiansMapper
const AbsolutePVCoordinates = @jimport org.orekit.utils.AbsolutePVCoordinates
const AbsolutePartialDerivativesEquations = @jimport org.orekit.propagation.numerical.AbsolutePartialDerivativesEquations
const AbstractGradientConverter = @jimport org.orekit.propagation.integration.AbstractGradientConverter
const AbstractIntegratedPropagator = @jimport org.orekit.propagation.integration.AbstractIntegratedPropagator
const AbstractJacobiansMapper = @jimport org.orekit.propagation.integration.AbstractJacobiansMapper
const AbstractPropagator = @jimport org.orekit.propagation.AbstractPropagator
const AdditionalEquations = @jimport org.orekit.propagation.integration.AdditionalEquations
const AdditionalStateProvider = @jimport org.orekit.propagation.AdditionalStateProvider
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BoundedPropagator = @jimport org.orekit.propagation.BoundedPropagator
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const EpochDerivativesEquations = @jimport org.orekit.propagation.numerical.EpochDerivativesEquations
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldAbstractIntegratedPropagator = @jimport org.orekit.propagation.integration.FieldAbstractIntegratedPropagator
const FieldAbstractPropagator = @jimport org.orekit.propagation.FieldAbstractPropagator
const FieldAdditionalEquations = @jimport org.orekit.propagation.integration.FieldAdditionalEquations
const FieldAdditionalStateProvider = @jimport org.orekit.propagation.FieldAdditionalStateProvider
const FieldBoundedPropagator = @jimport org.orekit.propagation.FieldBoundedPropagator
const FieldEventDetector = @jimport org.orekit.propagation.events.FieldEventDetector
const FieldNumericalPropagator = @jimport org.orekit.propagation.numerical.FieldNumericalPropagator
const FieldODEIntegrator = @jimport org.hipparchus.ode.FieldODEIntegrator
const FieldOrbit = @jimport org.orekit.orbits.FieldOrbit
const FieldOrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitFixedStepHandler
const FieldOrekitStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitStepHandler
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeDerivativesEquations = @jimport org.orekit.propagation.numerical.FieldTimeDerivativesEquations
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const GLONASSNumericalPropagator = @jimport org.orekit.propagation.numerical.GLONASSNumericalPropagator
const GLONASSNumericalPropagator_Builder = @jimport org.orekit.propagation.numerical.GLONASSNumericalPropagator$Builder
const GLONASSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GLONASSOrbitalElements
const Gradient = @jimport org.hipparchus.analysis.differentiation.Gradient
const JString = @jimport java.lang.String
const JacobiansMapper = @jimport org.orekit.propagation.numerical.JacobiansMapper
const List = @jimport java.util.List
const NumericalGradientConverter = @jimport org.orekit.propagation.numerical.NumericalGradientConverter
const NumericalPropagator = @jimport org.orekit.propagation.numerical.NumericalPropagator
const ODEIntegrator = @jimport org.hipparchus.ode.ODEIntegrator
const Object = @jimport java.lang.Object
const Orbit = @jimport org.orekit.orbits.Orbit
const OrbitType = @jimport org.orekit.orbits.OrbitType
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const PVCoordinates = @jimport org.orekit.utils.PVCoordinates
const ParameterDriversList = @jimport org.orekit.utils.ParameterDriversList
const PartialDerivativesEquations = @jimport org.orekit.propagation.numerical.PartialDerivativesEquations
const PositionAngle = @jimport org.orekit.orbits.PositionAngle
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeDerivativesEquations = @jimport org.orekit.propagation.numerical.TimeDerivativesEquations
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export AbsoluteJacobiansMapper
export AbsolutePartialDerivativesEquations
export EpochDerivativesEquations
export FieldNumericalPropagator
export FieldTimeDerivativesEquations
export GLONASSNumericalPropagator
export JacobiansMapper
export NumericalGradientConverter
export NumericalPropagator
export PartialDerivativesEquations
export TimeDerivativesEquations
export add_additional_equations
export add_additional_state_provider
export add_event_detector
export add_force_model
export add_kepler_contribution
export add_mass_derivative
export add_non_keplerian_acceleration
export clear_events_detectors
export compute_derivatives
export equals
export get_additional_state_dimension
export get_additional_state_providers
export get_all_force_models
export get_attitude_provider
export get_basic_dimension
export get_calls
export get_class
export get_events_detectors
export get_field
export get_frame
export get_free_state_parameters
export get_generated_ephemeris
export get_glonass_orbital_elements
export get_initial_state
export get_managed_additional_states
export get_mapper
export get_mode
export get_mu
export get_name
export get_orbit_type
export get_parameters
export get_parameters_jacobian
export get_position_angle_type
export get_pv_coordinates
export get_pv_in_pz90
export get_selected_parameters
export get_state
export get_state_jacobian
export hash_code
export init
export is_additional_state_managed
export is_initialize
export notify
export notify_all
export propagate
export remove_force_models
export reset_initial_state
export set_attitude_provider
export set_ephemeris_mode
export set_ignore_central_attraction
export set_initial_jacobians
export set_initial_state
export set_master_mode
export set_mu
export set_orbit_type
export set_position_angle_type
export set_reset_at_end
export set_slave_mode
export to_string
export tolerances
export wait

include("Cr3bpWrapper/Cr3bpWrapper.jl")
include("absolute_jacobians_mapper.jl")
include("absolute_partial_derivatives_equations.jl")
include("epoch_derivatives_equations.jl")
include("field_numerical_propagator.jl")
include("field_time_derivatives_equations.jl")
include("glonass_numerical_propagator.jl")
include("jacobians_mapper.jl")
include("numerical_gradient_converter.jl")
include("numerical_propagator.jl")
include("partial_derivatives_equations.jl")
include("time_derivatives_equations.jl")

end


text:
"""
    InnerProduct{F,S}

A struct representing a linear form on a Hilbert space as an inner product ``⟨w, x⟩``.
"""
struct InnerProduct{F,S}
    w::F
    solver::S
end

InnerProduct(w) = InnerProduct(w, HCubatureJL())

function (ip::InnerProduct)(x)
    product_integral(ip.w, x, ip.solver)
end

Flux.@functor InnerProduct

function w_pullback(ip::InnerProduct{<:ParametricFun,S}, x, δy) where {S}
    ∇w(u) = to_grad(f -> f(u), ip.w)
    s = product_integral(∇w, x)
    return from_grad(ip.w, s * δy)
end

w_pullback(ip, x, δy) = ZeroTangent()

function ChainRulesCore.rrule(ip::InnerProduct, x)
    y = ip(x)
    function ip_pullback(ȳ)
        δw = w_pullback(ip, x, ȳ)
        īp = Tangent{InnerProduct}(; w= δw, solver=NoTangent())
        x̄ = t -> ȳ * ip.w(t)
        return īp, x̄
    end
    return y, ip_pullback
end


text:
using Test
using CollegeEntry, StructLH

ce = CollegeEntry;

function common_test(a :: T1) where T1 <: AbstractAdmissionsRule
    @testset "$a" begin
        # println("\n------------------");
        # println(a);
        nc = n_colleges(a);
        af = ce.make_test_admprob_fct_logistic(nc);

        # Want these increasing
        hsGpaPctV = 0.01 : 0.2 : 0.99;
        J = length(hsGpaPctV);

        @test validate_admissions(a)
        # println(a);
        @test StructLH.describe(a) isa Matrix{String};
        nSets = n_college_sets(a);
        @test nSets >= 1
        @test n_colleges(a) > 1

        prob_jsM = zeros(J, nSets);
        for (iSet, cs) in enumerate(a)
            @test isa(collect(cs), Vector{<: Integer})
            @test isequal(cs,  college_set(a, iSet))
            setProbV = prob_coll_set(a, af, iSet, hsGpaPctV);
            @test all(setProbV .>= 0.0)
            @test all(setProbV .<= 1.0)
            # Cannot compare size b/c the size of a scalar is `nothing`
            @test length(setProbV) == length(hsGpaPctV)
            prob_jsM[:, iSet] = setProbV;

            # for j = 1 : J
            #     prob = CollegeStrat.prob_coll_set(a, iSet, hsGpaPctV[j]);
            #     @test isapprox(setProbV[j], prob, atol = 1e-6)
            # end
        end
        @test all(isapprox.(sum(prob_jsM, dims = 2), 1.0, atol = 1e-6))

        prob_jcM = admission_probs(a, af, hsGpaPctV);
        @test all(prob_jcM .>= 0.0)  &&  all(prob_jcM .<= 1.0)
        # Higher GPA students should be more likely to get into each college
        @test all(diff(prob_jcM, dims = 1) .>= -1e-8)
        # Better colleges should be harder to get into
        @test all(diff(prob_jcM, dims = 2) .<= 1e-8)
    end
end

function gpa_cutoff_test(a :: AdmissionsCutoff)
    @testset "GPA cutoff" begin
        hsGpaPctV = 0.01 : 0.3 : 0.99;
        highV = CollegeEntry.highest_college(a, hsGpaPctV);
        @test isa(highV, Vector{<: Integer})
        @test all(diff(highV) .>= 0)
        @test highV[1] == 1
        @test highV[end] == n_colleges(a);

        # Better students should have higher prob of "better" college set
        nSets = n_college_sets(a);
        admProbFct = CollegeEntry.make_test_admprob_fct_open(n_colleges(a));
        probV = prob_coll_set(a, admProbFct, nSets, hsGpaPctV);
        @test all(diff(probV) .>= 0)
        @test probV[end] > probV[1]
    end
end

@testset "Admission Rules" begin
    nc = 4;
    for a in [
        CollegeEntry.make_test_admissions_open(nc),
        CollegeEntry.make_test_admissions_cutoff(nc),
        CollegeEntry.make_test_admissions_onevar(nc)
        ]

        common_test(a);

        if isa(a, AdmissionsCutoff)
            gpa_cutoff_test(a);
        end
    end
end

# --------------


text:
export TR_Cub
function TR_Cub(h :: AbstractNLPModel;
               t₀ :: Float64 = h.meta.x0[1],
               tol :: Float64 = 1e-7,
               maxiter :: Int = 50,
               verbose :: Bool = true,
               eps1 ::Float64 = 0.2,
               eps2 :: Float64 = 0.8,
               red :: Float64 = 0.5,
               aug :: Float64 = 2.0,
               Δ :: Float64 = 1.0)

    (length(h.meta.x0) > 1) && warn("Not a 1-D problem ")
    t = t₀; iter = 0;
    fₖ = obj(h, [t])[1]
    gₖ = grad(h, [t])[1]

    secₖ = 1.0
    dN = -gₖ #pour la première itération
    A = 0.5; B = 0.0

    q(d) = fₖ + gₖ * d + A * d^2 + B * d^3

    verbose && @printf(" iter  t         gₖ          Δ        pred         ared\n")
    verbose && @printf(" %4d %7.2e  %7.2e  %7.2e \n", iter,t,gₖ,Δ)

    while ((abs(gₖ) > tol) & (iter < maxiter)) | (iter == 0)

        if q(Δ) < q(-Δ)
            d = Δ
        else
            d = -Δ
        end

        cub(t) = fₖ + gₖ * t + A * t^2 + B * t^3
        #dcub(t)=gₖ+2*A*t+3*B*t^2
        dR = roots(Poly([gₖ, 2 * A, 3 * B]))

        if isreal(dR)
          dR = real(dR)
          dN2 = dR[1]
          if length(dR) > 1
            if (cub(dR[2]) < cub(dR[1]))
              dN2 = dR[2]
            end
          end
        else
          dN2 = -gₖ * s / y
        end

        if (abs(dN2) < Δ) & (q(d) > q(dN2))
            d = dN2
        end

        if d == 0.0
            return (t, fₖ, norm(gₖ, Inf), iter + 1, false, false, :FailedDirectionComputation, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)
            #L'OUTIL DE CALCUL DES RACINES DU PACKAGE Polynomials ARRONDIE À 0 CERTAINES RACINES D'OÙ LA NOUVELLE CONDITION
        end

        ftestTR = obj(h, [t + d])[1]
        gtestTR = grad(h, [t + d])[1]

        pred = gₖ * d + A * d^2 + B * d^3

        if pred > - 1e-10
            ared = (gₖ + gtestTR) * d / 2
        else
            ared = ftestTR - fₖ
        end

        ratio = ared / pred
        if (ratio < eps1) & (abs(d) == Δ)
            Δ = red * Δ
        else
          #two points memory
          tpred = t; gkm1 = gₖ; fkm1 = fₖ
          t = t + d; gₖ = gtestTR; fₖ = ftestTR

          #cubic two points interpolation
          s = t - tpred; y = gₖ - gkm1

          α = -s
          z = gₖ + gkm1 + 3 * (fₖ - fkm1) / α
          discr = z^2 - gₖ * gkm1
          denom = gₖ + gkm1 + 2 * z
          B = 1 / 3 * (gₖ + gkm1 +2 * z) / (α * α)
          A = -(gₖ + z) / α

          if ratio > eps2
            Δ = aug * Δ
          end
        end

        iter += 1
        verbose && @printf(" %4d %7.2e  %7.2e  %7.2e %7.2e %7.2e\n", iter,t,gₖ,Δ,pred,ared)
    end

    status = :NotSolved
    (abs(gₖ) < tol) && (status = :Optimal)
    (iter >= maxiter) && (status = :Tired)
    tired = iter > maxiter
    optimal = abs(gₖ) < tol
    return (t, fₖ, norm(gₖ, Inf), iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)

end


text:
function DeSerializeAllCoord!(sys::System,array::Array{Float64,1})
    for i in 0:div(size(array)[1],3)-1
        sys.atoms.coords_all[i+1][1]=array[i*3+1]
        sys.atoms.coords_all[i+1][2]=array[i*3+2]
        sys.atoms.coords_all[i+1][3]=array[i*3+3]
    end
end

function DeSerializeAllVels!(sys::System,array::Array{Float64,1})
    for i in 0:div(size(array)[1],3)-1
        sys.atoms.vels_all[i+1][1]=array[i*3+1]
        sys.atoms.vels_all[i+1][2]=array[i*3+2]
        sys.atoms.vels_all[i+1][3]=array[i*3+3]
    end
end

function DeSerializeAllBonds!(sys::System,array::Array{Int64,1})
    for i in 0:div(size(array)[1],3)-1
        sys.atoms.bonds_all[i+1][1]=array[i*3+1]
        sys.atoms.bonds_all[i+1][2]=array[i*3+2]
        sys.atoms.bonds_all[i+1][3]=array[i*3+3]
    end
end

function DeSerializeAllMasses!(sys::System,array::Array{Float64,1})
    for i in 0:div(size(array)[1],1)-1
        sys.atoms.masses_all[i+1]=array[i+1]
    end
end

function DeSerializeAllTypes!(sys::System,array::Array{Int64,1})
    for i in 0:div(size(array)[1],1)-1
        sys.atoms.types_all[i+1]=array[i+1]
    end
end

function DeInitMolecules(sys::System)
    molecules=Dict()
    index_atom=1
    index_bond=1
    for comm_mole_i in sys.mpi.mole_index_comms
        for mole_i in comm_mole_i

            molecules[string(mole_i)]=Dict()

            molecules[string(mole_i)]["bonds"]=Dict()
            bond_mole_length=sys.mpi.bond_moles_length[mole_i]

            for bondi in 1:bond_mole_length
                bond_=zeros(Int64,3)
                bond_[1]=sys.atoms.bonds_all[index_bond][1]
                bond_[2]=sys.atoms.bonds_all[index_bond][2] - index_atom + 1
                bond_[3]=sys.atoms.bonds_all[index_bond][3] - index_atom + 1
                molecules[string(mole_i)]["bonds"][string(bondi)]=bond_
                index_bond=index_bond+1
            end

            molecules[string(mole_i)]["atoms"]=Dict()
            mole_length=sys.mpi.moles_length[mole_i]
            for atomi in 1:mole_length
                molecules[string(mole_i)]["atoms"][string(atomi)]=Dict( "type"=>sys.atoms.types_all[index_atom],
                                                                        "mass"=>sys.atoms.masses_all[index_atom],
                                                                        "coords"=>sys.atoms.coords_all[index_atom],
                                                                        "vels"=>sys.atoms.vels_all[index_atom])
                index_atom=index_atom+1
            end

        end
    end
    return molecules
end


text:
function encode_unit(output, counter, last_char)
    if last_char != '1'
        if counter == 1
            output *= last_char
        else
            output *= string(counter) * last_char
        end
    end
    return output
end

function encode(s)

    output = ""
    counter = 1
    last_char = '1'

    for c in s
        if c != last_char
            output = encode_unit(output, counter, last_char)
            counter = 1
            last_char = c
        else
            counter += 1
        end
    end

    output = encode_unit(output, counter, last_char)

    return output
end

function decode(s)
    output = ""
    substr = ""
    for c in s
        substr *= c
        if !isdigit(c)
            if length(substr) == 1
                output *= substr
                substr = ""
            else
                target = last(substr)
                num = Meta.parse(chop(substr, tail=1))
                output *= target ^ num
            end
            substr = ""
        end
    end

    return output
end

decode("2A3B4C")

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
module StrideArrays

# Write your package code here.

using VectorizationBase, ArrayInterface,
    SLEEFPirates, VectorizedRNG,
    LoopVectorization, LinearAlgebra,
    Random#, StackPointers#,
    # SpecialFunctions # Perhaps there is a better way to support erf?

using VectorizationBase: align, AbstractStridedPointer, AbstractSIMDVector, StridedPointer, gesp, pause, zstridedpointer, val_dense_dims, preserve_buffer
using LoopVectorization: CloseOpen
using Static: StaticInt, Zero, One, StaticBool, True, False
using ArrayInterface: size, strides, offsets, indices,
    static_length, axes, dense_dims, stride_rank, StrideIndex
using StrideArraysCore: AbstractStrideArray, AbstractStrideMatrix, AbstractStrideVector,
    AbstractPtrStrideArray, PtrArray, static_expr, rank_to_sortperm,
    StrideArray, StrideVector, StrideMatrix, similar_layout,
    @gc_preserve

using Octavian
using Octavian: MemoryBuffer

export @StrideArray, @gc_preserve, # @Constant,
    AbstractStrideArray, AbstractStrideVector, AbstractStrideMatrix,
    StrideArray, StrideVector, StrideMatrix,
    PtrArray,# PtrVector, PtrMatrix,
    # ConstantArray, ConstantVector, ConstantMatrix, allocarray,
    matmul!, matmul_serial!, mul!, *ˡ, StaticInt,
    matmul, matmul_serial
# LazyMap,

include("rand.jl")
include("blas.jl")
include("broadcast.jl")
include("miscellaneous.jl")

# include("precompile.jl")
# _precompile_()

end


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

JLLWrappers.@generate_main_file_header("Cuba")
JLLWrappers.@generate_main_file("Cuba", UUID("3bed1096-5ab2-53a1-92e2-6c1cc31d0f4b"))
end  # module Cuba_jll


text:
import JuMP, Cbc, XLSX

# make sure u set the working directory to current files folder
data = XLSX.readxlsx("../../data/data3_6.xlsx")

#General
General = data["General"];
l = General["B1"]
num_providers = General["B2"]
num_clients = General["B3"]

#Providers
Providers = data["Providers"];
provCostString  = string("A2:A", num_providers+1)
provStartString = string("B2:B", num_providers+1)
f = Providers[provCostString]
w = Providers[provStartString]

#Providers
Clients = data["Clients"];
clientStartString  = string("A2:A", num_clients+1)
clientDurationString  = string("B2:B", num_clients+1)
s_times = Clients[clientStartString]
d = Clients[clientDurationString]

#Distances
Distances = data["Distances"];
t = Distances[Distances.dimension]

# need to convert to floats to use 'sortperm' below
w = convert(Array{Float64,1}, w[:,1])
s_times = convert(Array{Float64,1}, s_times[:,1])
d = convert(Array{Float64,1}, d[:,1])
f = convert(Array{Float64,1}, f[:,1])

client_sorted = sortperm(s_times)  # client start times sorted (in index)

arr = [[f[i], w[i]] for i in 1:length(w)]
provider_sorted = sortperm(f) # providers cost and start time sorted

# heuristic function that computes feasible routes (but likely not optimal, i.e. minimum cost)
function heuristics_result(provider_sorted, client_sorted, w, s_times, t, d)
    past_client = zeros(0)
    num_served = zeros(num_providers, 2)
    num_served[:,1] = provider_sorted

    for m in provider_sorted  # for each index in sorted provider, and then for each client
        current_location = 1  # 1 is at depot(14HS) currently
        for n in client_sorted
            # only proceed if this client has not been assigned to a provider yet.
            if n ∉ past_client
                currentTime_provider = w[m]
                currentTime_client = s_times[n]
                travelTime = t[current_location, n+1] # n+1 is nth client in the t matrix

                #time to travel from provider's current location to the client plus the time when the provider is free.
                time_to_arrive = travelTime + currentTime_provider

                #Time to serve client "n".
                duration = d[n]
                #Time finished serving client "n".
                finish_time = currentTime_client + duration

                #Ensure the provider can arrive in time (at or earlier than service start times)
                # Must return back the origin (14HS) by hour l at latest.
                if currentTime_client >= time_to_arrive && finish_time <= l - t[n+1, 1]
                        println("Provider: ", m, " is assigned to client: ", n)
                        num_served[findfirst(isequal(m), num_served[:,1]), 2] += 1
                        w[m] = finish_time # update the time of provider m
                        current_location = n+1  # n+1 is nth client in the t-matrix
                        append!(past_client, n) # record service provided for client n
                end
            end
        end
        # end the algorithm if we have served all clients. so O(N*<M)?
        if length(past_client) == length(client_sorted)
            break
        end
    end
    println("\n")
    return convert(Array{Int,2}, num_served), convert(Array{Int,1}, past_client)
end

function compute_cost(num_served)
    cost = 0
    for i in 1:length(num_served[:,2])
        if num_served[:,2][i] > 0
            cost += f[num_served[:,1][i]]
        end
    end
    return cost
end

function is_feasible(past_client)
    return length(past_client) == length(client_sorted)
end

function get_heuristics_path(num_served, past_client, k)
    path = zeros(1)
    k_ind = findfirst(isequal(k), num_served[:,1])
    if k_ind == 1
        for i in 1:num_served[1,2]
            append!(path, past_client[i])
        end
    end
    if k_ind != 1
        start_ind = sum(num_served[1:k_ind-1,2]) + 1
        for i in start_ind:start_ind+num_served[k_ind,2]-1
            append!(path, past_client[i])
        end
    end
    return convert(Array{Int,1}, append!(path, zeros(1))) # convert to integer array
end

num_served, past_client = heuristics_result(provider_sorted, client_sorted, copy(w), s_times, t, d)

num_served

past_client

is_feasible(past_client)
compute_cost(num_served)

# k is the provider# you want to get the path for
k=4
path = get_heuristics_path(num_served, past_client, k)
path += convert(Array{Int,1}, ones(length(get_heuristics_path(num_served, past_client, k))))
println(path)