text:
function compute_xc_yc(tip)
    tip_xsiz, tip_ysiz = size(tip)
    xc = round(Int, tip_xsiz/2, RoundNearestTiesUp) - 1
    yc = round(Int, tip_ysiz/2, RoundNearestTiesUp) - 1
    return xc, yc
end

function ireflect(surface)
    surface2 = - surface[end:-1:1, end:-1:1]
    return surface2
end

function icmap(image, surface, tip; thresh=0.1)
    xc, yc = compute_xc_yc(tip)
    surf_xsiz, surf_ysiz = size(surface)
    tip_xsiz, tip_ysiz = size(tip)
    r = zeros(Int64, size(surface))
    x = y = 0
    num_contact = 0
    for i = 1:surf_xsiz
        for j = 1:surf_ysiz
            num_contact = 0
            pxmin = max(xc-tip_xsiz, 1-i)
            pymin = max(yc-tip_ysiz, 1-j)
            pxmax = min(xc-1, surf_xsiz-i)
            pymax = min(yc-1, surf_ysiz-j)
            for px = pxmin:pxmax
                for py = pymin:pymax
                    temp1 = image[i, j] - tip[xc-px, yc-py]
                    temp2 = surface[i+px, j+py]
                    if abs(temp1 - temp2) < thresh
                        num_contact += 1
                        x = i+px
                        y = j+py
                    end
                end
            end
            if num_contact == 1
                r[x, y] = 1
            end
        end
    end
    return r
end

function idilation(surface, tip)
    xc, yc = compute_xc_yc(tip)
    surf_xsiz, surf_ysiz = size(surface)
    tip_xsiz, tip_ysiz = size(tip)
    r = similar(surface)
    for i = 1:surf_xsiz
        for j = 1:surf_ysiz
            pxmin = max(i-surf_xsiz, -xc+1)
            pymin = max(j-surf_ysiz, -yc+1)
            pxmax = min(i-1, -xc+tip_xsiz)
            pymax = min(j-1, -yc+tip_ysiz)
            dil_max = surface[i-pxmin, j-pymin] + tip[xc+pxmin, yc+pymin]
            for px = pxmin:pxmax
                for py = pymin:pymax
                    temp = surface[i-px, j-py] + tip[xc+px, yc+py]
                    dil_max = max(temp, dil_max)
                end
            end
            r[i, j] = dil_max
        end
    end
    return r
end

function ierosion(image, tip)
    xc, yc = compute_xc_yc(tip)
    im_xsiz, im_ysiz = size(image)
    tip_xsiz, tip_ysiz = size(tip)
    r = similar(image)
    for i = 1:im_xsiz
        for j = 1:im_ysiz
            pxmin = max(-i+1, -xc+1)
            pymin = max(-j+1, -yc+1)
            pxmax = min(-i+im_xsiz, -xc+tip_xsiz)
            pymax = min(-j+im_ysiz, -yc+tip_ysiz)
            eros_min = image[i+pxmin, j+pymin] - tip[xc+pxmin, yc+pymin]
            for px = pxmin:pxmax
                for py = pymin:pymax
                    temp = image[i+px, j+py] - tip[xc+px, yc+py]
                    eros_min = min(temp, eros_min)
                end
            end
            r[i, j] = eros_min
        end
    end
    return r
end

function iopen(image, tip)
    r = ierosion(image, tip)
    r = idilation(r, tip)
    return r
end

function itip_estimate_point!(tip0, image, ixp, jxp; thresh=0.0)
    xc, yc = compute_xc_yc(tip0)
    im_xsiz, im_ysiz = size(image)
    tip_xsiz, tip_ysiz = size(tip0)
    interior = (tip_xsiz <= (ixp-1) ) & (ixp <= (im_xsiz-tip_xsiz)) & (tip_ysiz <= (jxp-1)) & (jxp <= (im_ysiz-tip_ysiz))
    count = 0
    if interior
        for ix = 0:(tip_xsiz-1)
            for jx = 0:(tip_ysiz-1)
                if (ix+1) == xc & (jx+1) == yc
                    continue
                end
                dil = - eltype(image)(Inf)
                for id = 0:(tip_xsiz-1)
                    for jd = 0:(tip_ysiz-1)
                        if ((image[ixp+1, jxp+1] - image[ixp+xc-id, jxp+yc-jd]) > 0.0) | (xc == (id+1) & yc == (jd+1))
                            continue
                        end
                        temp = image[ixp+ix-id+1, jxp+jx-jd+1] - image[ixp+1, jxp+1] + thresh
                        dil = max(dil, temp)
                    end
                end
                if isinf(-dil)
                    continue
                end
                if dil < tip0[ix+1, jx+1]
                    tip0[ix+1, jx+1] = dil
                    count += 1
                end
            end
        end
    end
    return count
end

function itip_estimate_iter!(tip0, image::Matrix{T}; thresh=0.0) where {T <: Number}
    xc, yc = compute_xc_yc(tip0)
    im_xsiz, im_ysiz = size(image)
    tip_xsiz, tip_ysiz = size(tip0)
    count = 0
    open_image = iopen(image, tip0)
    for ixp = (tip_xsiz-xc-1):(im_xsiz-xc-1)
        for jxp = (tip_ysiz-yc-1):(im_ysiz-yc-1)
            if (image[ixp+1, jxp+1] - open_image[ixp+1, jxp+1]) > thresh
                c = itip_estimate_point!(tip0, image, ixp, jxp, thresh=thresh)
                if c > 0
                    count += 1
                end
            end
        end
    end
    return count
end

function itip_estimate!(tip0, images::Vector{Any}; thresh=0.0)
    nframe = length(images)
    for iframe = 1:nframe
        @printf "Processing %d-th image\n" iframe
        itip_estimate!(tip0, images[iframe], thresh=0.0)
        @printf "\n"
    end
    return
end

function itip_estimate!(tip0, image::Matrix{T}; thresh=0.0) where {T <: Number}
    iter = 0
    count = 1
    while count > 0
        iter += 1
        count = itip_estimate_iter!(tip0, image, thresh=thresh)
        @printf "Finished iteration %d\n" iter
        @printf "%d image locations produced refinement\n" count
    end
    return count
end

function idilation_pdiff(surface, tip)
    xc, yc = compute_xc_yc(tip)
    surf_xsiz, surf_ysiz = size(surface)
    tip_xsiz, tip_ysiz = size(tip)
    r = similar(surface)
    r_index = zeros(Int, size(r))
    for i = 1:surf_xsiz
        for j = 1:surf_ysiz
            pxmin = max(i-surf_xsiz, -xc+1)
            pymin = max(j-surf_ysiz, -yc+1)
            pxmax = min(i-1, -xc+tip_xsiz)
            pymax = min(j-1, -yc+tip_ysiz)
            dil_max = surface[i-pxmin, j-pymin] + tip[xc+pxmin, yc+pymin]
            r_index[i, j] = (yc+pymin-1)*size(tip, 1) + (xc+pxmin)
            for px = pxmin:pxmax
                for py = pymin:pymax
                    temp = surface[i-px, j-py] + tip[xc+px, yc+py]
                    if temp > dil_max
                        dil_max = temp
                        r_index[i, j] = (yc+py-1)*size(tip, 1) + (xc+px)
                    end
                end
            end
            r[i, j] = dil_max
        end
    end
    return r, r_index
end

function ierosion_pdiff(image, tip)
    xc, yc = compute_xc_yc(tip)
    im_xsiz, im_ysiz = size(image)
    tip_xsiz, tip_ysiz = size(tip)
    r = similar(image)
    r_index = zeros(Int, size(r))
    for i = 1:im_xsiz
        for j = 1:im_ysiz
            pxmin = max(-i+1, -xc+1)
            pymin = max(-j+1, -yc+1)
            pxmax = min(-i+im_xsiz, -xc+tip_xsiz)
            pymax = min(-j+im_ysiz, -yc+tip_ysiz)
            eros_min = image[i+pxmin, j+pymin] - tip[xc+pxmin, yc+pymin]
            r_index[i, j] = (yc+pymin-1)*size(tip, 1) + (xc+pxmin)
            for px = pxmin:pxmax
                for py = pymin:pymax
                    temp = image[i+px, j+py] - tip[xc+px, yc+py]
                    if temp < eros_min
                        eros_min = temp
                        r_index[i, j] = (yc+py-1)*size(tip, 1) + (xc+px)
                    end
                end
            end
            r[i, j] = eros_min
        end
    end
    return r, r_index
end

function itip_least_squares!(tip0::Matrix{T}, image::Matrix{T}, surface::Matrix{T}; rate=0.01, max_convergence=0.1) where {T <: Number}
    tip_xsiz, tip_ysiz = size(tip0)
    icount = 0
    check_convergence = T(Inf)
    tip0_old = deepcopy(tip0)
    while check_convergence > max_convergence
        xc, yc = MDToolbox.compute_xc_yc(tip0)
        r, r_index = MDToolbox.idilation_pdiff(surface, tip0)
        for ix = 1:tip_xsiz
            for iy = 1:tip_ysiz
                if (ix == xc) & (iy == yc)
                    continue
                end
                id = r_index .== ((iy-1)*size(tip0, 1) + ix)
                if any(id)
                    tip0[ix, iy] -= sum(r[id] .- image[id]) * rate
                end
            end
        end
        icount += 1
        check_convergence = maximum(abs.(tip0_old .- tip0))
        tip0_old .= tip0
        if mod(icount, 10) == 0
            #@printf "iteration %d : %f\n" icount sum((idilation(surface, tip0) .- image).^2)
            @printf "iteration %d : %f\n" icount check_convergence
        end
    end
end

#function itip_em!(tip0::Matrix{T}, images::Array{Matrix{Float64}, 1}; rate=0.01, max_convergence=0.1, nsample=1000, max_iteration=1000) where {T <: Number}
function itip_em!(tip0::Matrix{T}, images::Vector{Any}; rate=0.01, max_convergence=0.1, nsample=1000, max_iteration=1000) where {T <: Number}
    xc, yc = compute_xc_yc(tip0)
    tip_xsiz, tip_ysiz = size(tip0)
    nframe = length(images)
    icount = 0
    check_convergence = T(Inf)
    tip0_old = deepcopy(tip0)
    delta = similar(tip0)
    surfaces = deepcopy(images)
    surfaces_sample = deepcopy(images)

    while (check_convergence > max_convergence) & (icount < max_iteration)
        for iframe = 1:nframe
            surfaces[iframe] .= ierosion(images[iframe], tip0)
        end

        # E-step
        ids = []
        vecs = []
        for iframe = 1:nframe
            cmap = icmap(images[iframe], surfaces[iframe], tip0)
            id = cmap .< 0.5
            vec_uncertain = surfaces_sample[iframe][id]
            vec = []
            for isample = 1:nsample
                push!(vec, vec_uncertain .* rand(eltype(vec_uncertain), length(vec_uncertain)))
            end
            push!(ids, deepcopy(id))
            push!(vecs, deepcopy(vec))
        end

        # M-step
        tip0_old .= tip0
        for istep = 1:100
            delta .= 0.0
            for iframe = 1:nframe
                for isample = 1:nsample
                    surface = surfaces[iframe]
                    surface[ids[iframe]] .= vecs[iframe][isample]
                    r, r_index = MDToolbox.idilation_pdiff(surface, tip0)
                    for ix = 1:tip_xsiz
                        for iy = 1:tip_ysiz
                            if (ix == xc) & (iy == yc)
                                continue
                            end
                            id = r_index .== ((iy-1)*size(tip0, 1) + ix)
                            if any(id)
                                delta[ix, iy] -= sum(r[id] .- images[iframe][id]) * rate / nsample
                            end
                        end
                    end
                end
            end
            tip0 .+= delta
            tip0 .= min.(tip0, 0.0)
        end
        icount += 1
        check_convergence = maximum(abs.(tip0_old .- tip0))
        if mod(icount, 1) == 0
            #@printf "iteration %d : %f\n" icount sum((idilation(surface, tip0) .- image).^2)
            @printf "iteration %d : %f\n" icount check_convergence
        end
    end
    return delta
end

function itip_least_squares!(tip0, image::Matrix{T}; thresh=0.1, rate=0.1) where {T <: Number}
    images = []
    push!(images, image)
    itip_least_squares!(tip0, images, thresh=thresh, rate=rate)
end

function itip_least_squares!(tip0, images::Vector{Any}; thresh=0.1, rate=0.1, nstep=100, lambda=0.0)
    #@show sum((idilation(surf, tip0) .- image).^2)
    xc, yc = MDToolbox.compute_xc_yc(tip0)
    tip0 .-= 10^(-8).*rand(Float64, size(tip0))
    tip0[xc, yc] = 0.0
    nframe = ceil(Int, length(images)*0.8)
    mframe = length(images)
    d = zeros(eltype(tip0), size(tip0))
    loss_train_array = []
    loss_val_array = []
    for i = 1:nstep
        d .= eltype(tip0)(0.0)
        #Threads.@threads for iframe = 1:nframe
        for iframe = 1:nframe
            s, e_index = ierosion_pdiff(images[iframe], tip0)
            r, d_index = idilation_pdiff(s, tip0)
            for ix = 1:size(tip0, 1)
                for iy = 1:size(tip0, 2)
                    if (ix == xc) & (iy == yc)
                        continue
                    end
                    id_e = e_index .== ((iy-1)*size(tip0, 1) + ix)
                    id_d = d_index .== ((iy-1)*size(tip0, 1) + ix)
                    #if mod(i, 2) == 0
                    #    d[ix, iy] += - sum((r .- images[iframe]) .* id_d) .* rate
                    #else
                    #    d[ix, iy] += - sum((r .- images[iframe]) .* (- id_e)) .* rate
                    #end
                    d[ix, iy] += - sum((r .- images[iframe]) .* (id_d .- id_e)) .* rate
                    #d[ix, iy] += - sum((r .- images[iframe]) .* (+id_d) .* (-id_e)) .* rate
                    #d[ix, iy] += lambda * rate
                end
            end
        end
        #for ix = 1:size(tip0, 1)
        #    for iy = 1:size(tip0, 2)
        #        if abs(d[ix, iy]) > 10.0
        #            d[ix, iy] = 10.0*d[ix, iy] / abs(d[ix, iy])
        #        end
        #    end
        #end
        #d .= min.(d, 0.0)
        tip0 .+= d
        tip0 .= min.(tip0, 0.0)
        #tip0 .= max.(tip0, -100.0)
        #tip0 .= tip0 .- tip0[xc, yc]

        #tip0_old = tip0
        #for ix = 2:(size(tip0, 1)-1)
        #    for iy = 2:(size(tip0, 2)-1)
        #        tip0[ix, iy] = 0.25 * (tip0_old[ix, iy-1] + tip0_old[ix-1, iy] + tip0_old[ix+1, iy] + tip0_old[ix, iy+1])
        #    end
        #end
        #tip0 .= tip0 .- tip0[xc, yc]

        if mod(i, 10) == 0
            loss_train = 0.0
            for iframe = 1:nframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                loss_train += sum((r .- images[iframe]).^2)
            end
            loss_train += lambda * sum(tip0.^2)
            #loss_train /= nframe

            loss_val = 0.0
            for iframe = (nframe+1):mframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                loss_val += sum((r .- images[iframe]).^2)
            end
            if nframe != mframe
                loss_val /= (mframe - nframe)
            end

            println("step $(i): loss_train = $(loss_train)   loss_validation = $(loss_val)")
            push!(loss_train_array, loss_train)
            push!(loss_val_array, loss_val)
        end
    end
    return loss_train_array, loss_val_array
end

function itip_least_squares2!(tip0, images::Vector{Any}; thresh=0.1, rate=0.1, nstep=100)
    #@show sum((idilation(surf, tip0) .- image).^2)
    xc, yc = MDToolbox.compute_xc_yc(tip0)
    tip0 .-= 10^(-8).*rand(Float64, size(tip0))
    tip0[xc, yc] = 0.0
    #tip0[end, end] = -100.0
    nframe = ceil(Int, length(images)*0.8)
    mframe = length(images)
    d = zeros(eltype(tip0), size(tip0))
    loss_train_array = []
    loss_val_array = []
    for i = 1:nstep
        d .= eltype(tip0)(0.0)
        for iframe = 1:nframe
            s, e_index = ierosion_pdiff(images[iframe], tip0)
            #r, d_index = idilation_pdiff(s, tip0)
            s2, e2_index = idilation_pdiff(s, tip0)
            r2, d2_index = ierosion_pdiff(s2, tip0)
            @show sum(abs.(s .- r2))
            return
            for ix = 1:size(tip0, 1)
                for iy = 1:size(tip0, 2)
                    if (ix == xc) & (iy == yc)
                        continue
                    end
                    id_e = e_index .== ((iy-1)*size(tip0, 1) + ix)
                    #id_d = d_index .== ((iy-1)*size(tip0, 1) + ix)
                    id_e2 = e2_index .== ((iy-1)*size(tip0, 1) + ix)
                    id_d2 = d2_index .== ((iy-1)*size(tip0, 1) + ix)
                    #if mod(i, 2) == 0
                    #    d[ix, iy] += - sum((r .- images[iframe]) .* id_d) .* rate
                    #else
                    #    d[ix, iy] += - sum((r .- images[iframe]) .* (- id_e)) .* rate
                    #end
                    #d[ix, iy] += - sum((r .- images[iframe]) .* (id_d .- id_e)) .* rate
                    d[ix, iy] += - sum((r2 .- s) .* (id_e2 .- id_d2)) .* rate
                    #d[ix, iy] += 10.0*(d[xc, yc] - d[ix, iy]) * rate
                end
            end
        end
        #for ix = 1:size(tip0, 1)
        #    for iy = 1:size(tip0, 2)
        #        if abs(d[ix, iy]) > 10.0
        #            d[ix, iy] = 10.0*d[ix, iy] / abs(d[ix, iy])
        #        end
        #    end
        #end
        #d .= min.(d, 0.0)
        tip0 .+= d
        tip0 .= min.(tip0, 0.0)
        #tip0 .= tip0 .- tip0[xc, yc]

        #tip0_old = tip0
        #for ix = 2:(size(tip0, 1)-1)
        #    for iy = 2:(size(tip0, 2)-1)
        #        tip0[ix, iy] = 0.25 * (tip0_old[ix, iy-1] + tip0_old[ix-1, iy] + tip0_old[ix+1, iy] + tip0_old[ix, iy+1])
        #    end
        #end
        #tip0 .= tip0 .- tip0[xc, yc]

        if mod(i, 10) == 0
            loss_train = 0.0
            for iframe = 1:nframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                r2, d_index = ierosion_pdiff(r, tip0)
                loss_train += sum((r2 .- s).^2)
            end
            loss_train /= nframe

            loss_val = 0.0
            for iframe = (nframe+1):mframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                loss_val += sum((r .- images[iframe]).^2)
            end
            if nframe != mframe
                loss_val /= (mframe - nframe)
            end

            println("step $(i): loss_train = $(loss_train)   loss_validation = $(loss_val)")
            push!(loss_train_array, loss_train)
            push!(loss_val_array, loss_val)
        end
    end
    return loss_train_array, loss_val_array
end

function itip_least_squares_adam2!(tip0, images::Vector{Any}; thresh=0.1, learning_rate=0.001)
    xc, yc = MDToolbox.compute_xc_yc(tip0)
    tip0 .-= 10^(-8).*rand(Float64, size(tip0))
    tip0[xc, yc] = 0.0
    nframe = length(images)
    BETA1 = 0.9
    BETA2 = 0.999
    EPS = 1e-08
    mt = similar(tip0)
    mt .= 0.0
    vt = similar(tip0)
    vt .= 0.0
    delta = similar(tip0)
    delta .= 0.0
    for i = 1:1000
        #lambda = 0.1
        #ss = tip0
        delta .= 0.0
        for iframe = 1:nframe
        #for iframe in randperm(MersenneTwister(1234), nframe)
            s1, e1_index = ierosion_pdiff(images[iframe], tip0)
            r1, d1_index = idilation_pdiff(s1, tip0)
            s2, e2_index = ierosion_pdiff(r1, tip0)
            r2, d2_index = idilation_pdiff(s2, tip0)
            s3, e3_index = ierosion_pdiff(r2, tip0)
            r3, d3_index = idilation_pdiff(s3, tip0)
            for ix = 1:size(tip0, 1)
                for iy = 1:size(tip0, 2)
                    if (ix == xc) & (iy == yc)
                        continue
                    end
                    id1_e = e1_index .== ((iy-1)*size(tip0, 1) + ix)
                    id1_d = d1_index .== ((iy-1)*size(tip0, 1) + ix)
                    id2_e = e2_index .== ((iy-1)*size(tip0, 1) + ix)
                    id2_d = d2_index .== ((iy-1)*size(tip0, 1) + ix)
                    id3_e = e3_index .== ((iy-1)*size(tip0, 1) + ix)
                    id3_d = d3_index .== ((iy-1)*size(tip0, 1) + ix)
                    delta[ix, iy] += sum((r3 .- images[iframe]) .* (id1_d .- id1_e .+ id2_d .- id2_e .+ id3_d .- id3_e)) .* learning_rate
                    #if (ix >= 2) & (ix <= (size(tip0, 1)-1))
                    #    if (iy >= 2) & (iy <= (size(tip0, 2)-1))
                    #        delta[ix, iy] -= lambda*((ss[ix, iy] - ss[ix, iy+1]) + (ss[ix, iy] - ss[ix, iy-1]) + (ss[ix, iy] - ss[ix+1, iy]) + (ss[ix, iy] - ss[ix-1, iy]))
                    #    end
                    #end
                end
            end
        end
        @. mt = BETA1 * mt + (1 - BETA1) * delta
        @. vt = BETA2 * vt + (1 - BETA2) * delta^2
        @. delta = mt / (1.0 - BETA1^i) / (sqrt(vt / (1.0 - BETA2^i)) + EPS) * learning_rate
        #delta .= max.(delta, 0.0)
        tip0 .= tip0 .- delta
        tip0 .= min.(tip0, 0.0)
        if mod(i, 100) == 0
            loss = 0.0
            for iframe = 1:nframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                loss += sum((r .- images[iframe]).^2)
            end
            println("step $(i): loss = $(loss)")
        end
    end
    return
end

function itip_least_squares_adam!(tip0, images::Vector{Any}; thresh=0.1, learning_rate=0.001)
    xc, yc = MDToolbox.compute_xc_yc(tip0)
    tip0 .-= 10^(-8).*rand(Float64, size(tip0))
    tip0[xc, yc] = 0.0
    nframe = length(images)
    BETA1 = 0.9
    BETA2 = 0.999
    EPS = 1e-08
    mt = similar(tip0)
    mt .= 0.0
    vt = similar(tip0)
    vt .= 0.0
    delta = similar(tip0)
    delta .= 0.0
    loss_array = []
    for i = 1:1000
        #lambda = 0.1
        #ss = tip0
        delta .= 0.0
        for iframe = 1:nframe
        #for iframe in randperm(MersenneTwister(1234), nframe)
            s, e_index = ierosion_pdiff(images[iframe], tip0)
            r, d_index = idilation_pdiff(s, tip0)
            for ix = 1:size(tip0, 1)
                for iy = 1:size(tip0, 2)
                    if (ix == xc) & (iy == yc)
                        continue
                    end
                    id_e = e_index .== ((iy-1)*size(tip0, 1) + ix)
                    id_d = d_index .== ((iy-1)*size(tip0, 1) + ix)
                    delta[ix, iy] += sum((r .- images[iframe]) .* (id_d .- id_e)) .* learning_rate
                    #if (ix >= 2) & (ix <= (size(tip0, 1)-1))
                    #    if (iy >= 2) & (iy <= (size(tip0, 2)-1))
                    #        delta[ix, iy] -= lambda*((ss[ix, iy] - ss[ix, iy+1]) + (ss[ix, iy] - ss[ix, iy-1]) + (ss[ix, iy] - ss[ix+1, iy]) + (ss[ix, iy] - ss[ix-1, iy]))
                    #    end
                    #end
                end
            end
        end
        @. mt = BETA1 * mt + (1 - BETA1) * delta
        @. vt = BETA2 * vt + (1 - BETA2) * delta^2
        @. delta = mt / (1.0 - BETA1^i) / (sqrt(vt / (1.0 - BETA2^i)) + EPS) * learning_rate
        #delta .= max.(delta, 0.0)
        tip0 .= tip0 .- delta
        #tip0 .-= 10^(-8).*rand(Float64, size(tip0))
        tip0 .= min.(tip0, 0.0)
        #tip0 .= tip0 .- tip0[xc, yc]
        if mod(i, 10) == 0
            loss = 0.0
            for iframe = 1:nframe
                s, e_index = ierosion_pdiff(images[iframe], tip0)
                r, d_index = idilation_pdiff(s, tip0)
                loss += sum((r .- images[iframe]).^2)
            end
            println("step $(i): loss = $(loss)")
            push!(loss_array, loss)
        end
    end
    return loss_array
end

"""
afmize

simulate afm image from a given structure

example::
pdb = readpbd("1ake.pdb");
decenter!(pdb);
afm_image = afmize(pdb, (150.0, 150.0), (16, 16));
heatmap(afm_image) # using Plots
"""
function afmize_alpha(ta::TrjArray, box, npixel)::Matrix{Float64}
    rpixel = box ./ npixel
    nx = floor.(Int, (ta.xyz[1, 1:3:end] .+ 0.5*box[1]) ./ rpixel[1]);
    ny = floor.(Int, (ta.xyz[1, 2:3:end] .+ 0.5*box[2]) ./ rpixel[2]);
    #@show nx
    #@show ny
    afm_image = zeros(Float64, npixel)
    z_min = minimum(ta.xyz[1, 3:3:end])
    z = ta.xyz[1, 3:3:end] .- z_min
    for iatom in 1:ta.natom
        if (1 <= nx[iatom] <= npixel[1]) && (1 <= ny[iatom] <= npixel[2])
            if z[iatom] > afm_image[nx[iatom], ny[iatom]]
                afm_image[nx[iatom], ny[iatom]] = z[iatom]
            end
        end
    end
    return afm_image
end

function defaultParameters()
    return Dict(
     "H" => 1.20,
    "HE" => 1.40,
     "B" => 1.92,
     "C" => 1.70,
     "CA" => 1.70,
     "CB" => 1.70,
     "CG" => 1.70,
     "CG1" => 1.70,
     "CG2" => 1.70,
     "CG3" => 1.70,
     "CD" => 1.70,
     "CD1" => 1.70,
     "CD2" => 1.70,
     "CD3" => 1.70,
     "CZ" => 1.70,
     "CZ1" => 1.70,
     "CZ2" => 1.70,
     "CZ3" => 1.70,
     "CD" => 1.70,
     "CD1" => 1.70,
     "CD2" => 1.70,
     "CD3" => 1.70,
     "CE" => 1.70,
     "CE1" => 1.70,
     "CE2" => 1.70,
     "CE3" => 1.70,
     "CH" => 1.70,
     "CH1" => 1.70,
     "CH2" => 1.70,
     "CH3" => 1.70,
     "N" => 1.55,
     "NE" => 1.55,
     "NZ" => 1.55,
     "ND1" => 1.55,
     "ND2" => 1.55,
     "NE1" => 1.55,
     "NE2" => 1.55,
     "NH1" => 1.55,
     "NH2" => 1.55,
      "O" => 1.52,
     "OH" => 1.52,
     "OG" => 1.52,
     "OE1" => 1.52,
     "OE2" => 1.52,
     "OG1" => 1.52,
     "OG2" => 1.52,
     "OD1" => 1.52,
     "OD2" => 1.52,
     "OXT" => 1.52,
     "F" => 1.47,
    #"NE" => 1.54,
    "MG" => 1.73,
    "AL" => 1.84,
    "SI" => 2.10,
     "P" => 1.80,
     "S" => 1.80,
     "SD" => 1.80,
     "SG" => 1.80,
    "CL" => 1.75,
    "AR" => 1.88,
     "K" => 2.75,
    "CA" => 2.31,
    "CYS" => 2.75,
    "PHE" => 3.2,
    "LEU" => 3.1,
    "TRP" => 3.4,
    "VAL" => 2.95,
    "ILE" => 3.1,
    "MET" => 3.1,
    "HIS" => 3.05,
    "HSD" => 3.05,
    "TYR" => 3.25,
    "ALA" => 2.5,
    "GLY" => 2.25,
    "PRO" => 2.8,
    "ASN" => 2.85,
    "THR" => 2.8,
    "SER" => 2.6,
    "ARG" => 3.3,
    "GLN" => 3.0,
    "ASP" => 2.8,
    "LYS" => 3.2,
    "GLU" => 2.95,
    )
end

mutable struct Point2D{T<:Real}
    x::T
    y::T
end

mutable struct Point3D{T<:Real}
    x::T
    y::T
    z::T
end

struct Probe
    x::Float64
    y::Float64
    r::Float64
    angle::Float64
end

mutable struct Sphere
    x::Float64
    y::Float64
    z::Float64
    r::Float64
end

mutable struct AfmizeConfig
    probeAngle::Float64        # Radian
    probeRadius::Float64
    range_min::Point2D
    range_max::Point2D
    resolution::Point2D
    atomRadiusDict::Dict{String, Float64}
end

function afmize!(tip::Matrix, config::AfmizeConfig)
    tip_xsiz, tip_ysiz = size(tip)
    xc, yc = compute_xc_yc(tip)
    #xc += 1
    #yc += 1
    @show xc, yc
    for ix = 1:tip_xsiz
        for iy = 1:tip_ysiz
            x = config.resolution.x * abs(ix - xc)
            y = config.resolution.y * abs(iy - yc)
            d = sqrt(x^2 + y^2)
            if d <= config.probeRadius
                z = sqrt(config.probeRadius^2 - d^2)
            else
                #theta = (0.5*pi) - (0.5*config.probeAngle)
                theta = (0.5*pi) - (config.probeAngle)
                z = - tan(theta) * (d - config.probeRadius)
            end
            tip[ix, iy] = z
        end
    end
    tip .= tip .- maximum(tip)
    return 0
end

"""
the following codes are licensed under the MIT license (Copyright (c) 2018-2021 Tohru Niina), see LICENSE.md
"""

function AfmizeConfig(r::Float64)
    return AfmizeConfig(10.0 * (pi / 180),
                        r,
                        Point2D(-90, -90),
                        Point2D(90, 90),
                        Point2D(15, 15),
                        defaultParameters())
end

function defaultConfig()
    return AfmizeConfig(10.0 * (pi / 180),
                        20.0,
                        Point2D(-90, -90),
                        Point2D(90, 90),
                        Point2D(15, 15),
                        defaultParameters())
end

function checkConfig(tra::TrjArray ,config::AfmizeConfig)::Union{String, Nothing}
    if (config.range_max.x - config.range_min.x) % config.resolution.x != 0
        return "resolution.x must divide range"
    end
    if (config.range_max.y - config.range_min.y) % config.resolution.y != 0
        return "resolution.y must divide range"
    end
    if config.probeAngle == 0
        return "probeAngle must be positive"
    end

    for name in tra.atomname
        if !haskey(config.atomRadiusDict, name)
            return "config doesn't know atom name $(name)"
        end
    end

    return nothing
end

# 底を0にする
function moveBottom(atoms::Array{Sphere})
    bottom = atoms[1].z - atoms[1].r
    for atom in atoms
        bottom = min(bottom, atom.z - atom.r)
    end
    for atom in atoms
        atom.z -= bottom
    end
end

# 針を球として見たとき、どこで衝突するかの計算
function calcCollisionAsSphere(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    dr = probe.r + atom.r
    if 0 < dr^2 - distXY^2
        return atom.z + sqrt(dr^2 - distXY^2) - probe.r
    else
        return 0.0
    end
end

# 針を円形推台として見たとき、どこで衝突するかの計算
function calcCollisionAsCircularThrusters(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    collisionDist = probe.r + atom.r * cos(probe.angle)
    if collisionDist < distXY
        return atom.z + atom.r * sin(probe.angle) - (distXY - collisionDist) / tan(probe.angle) - probe.r
    elseif probe.r < distXY
        return atom.z + sqrt(atom.r^2 - (distXY - probe.r)^2) - probe.r
    else
        return 0.0
    end
end

# 衝突計算をする範囲を返す
function calRectangle(min_point::Point2D, max_point::Point2D, config::AfmizeConfig)
    if max_point.x <= config.range_min.x || max_point.y <= config.range_min.y
        return nothing
    end
    if config.range_max.x <= min_point.x || config.range_max.y <= min_point.y
        return nothing
    end

    min_point = Point2D(max(min_point.x, config.range_min.x),
                        max(min_point.y, config.range_min.y))
    max_point = Point2D(min(max_point.x, config.range_max.x),
                        min(max_point.y, config.range_max.y))
    width  = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    return (max(1, floor(Int, (min_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            max(1, floor(Int, (min_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)),
            min(width , floor(Int, (max_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            min(height, floor(Int, (max_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)))
end

function calRectangle_sphere(atom::Sphere, config::AfmizeConfig)
    min_point = Point2D(atom.x - atom.r - config.probeRadius,
                        atom.y - atom.r - config.probeRadius)
    max_point = Point2D(atom.x + atom.r + config.probeRadius,
                        atom.y + atom.r + config.probeRadius)

    return calRectangle(min_point, max_point, config)
end

function calRectangle_circularThrusters(atom::Sphere, config::AfmizeConfig)
    dist_xy = (tan(config.probeAngle) * (atom.z + atom.r * sin(config.probeAngle) - config.probeRadius)
                + config.probeRadius + atom.r * cos(config.probeAngle))
    min_point = Point2D(atom.x - dist_xy,
                        atom.y - dist_xy)
    max_point = Point2D(atom.x + dist_xy,
                        atom.y + dist_xy)

    return calRectangle(min_point, max_point, config)
end

function surfing(t::TrjArray, config::AfmizeConfig)
    message = checkConfig(t, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    stage = zeros(height, width)

    #radius_max = maximum(t.radius)
    radius = zeros(Float64, t.natom)
    for iatom = 1:t.natom
        radius[iatom] = config.atomRadiusDict[t.atomname[iatom]]
    end

    x = t.xyz[1, 1:3:end]
    y = t.xyz[1, 2:3:end]
    z = t.xyz[1, 3:3:end]
    z .= z .- minimum(z)
    for w in 1:width
        grid_x = config.range_min.x + (w-0.5) * config.resolution.x
        dx = abs.(grid_x .- x)
        index_x = dx .< radius
        for h in 1:height
            grid_y = config.range_min.y + (h-0.5) * config.resolution.y
            dy = abs.(grid_y .- y)
            index_y = dy .< radius
            index = index_x .& index_y
            if any(index)
                d = sqrt.(dx[index].^2 + dy[index].^2)
                r = radius[index]
                index2 = d .< r
                if any(index2)
                    z_surface = z[index][index2] .+ sqrt.(r[index2].^2 .- d[index2].^2)
                    stage[h, w] = maximum(z_surface)
                end
            end
        end
    end

    return stage
end

function afmize(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.xyz[1, 3*(i-1)+1], tra.xyz[1, 3*(i-1)+2], tra.xyz[1, 3*(i-1)+3],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)
    probes = [Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                    config.range_min.y + (h-0.5) * config.resolution.y,
                    config.probeRadius, config.probeAngle)
             for h in 1:height, w in 1:width]

    # 各原子事に計算する価値のある矩形を求めて、その中で探索を行う
    for atom in atoms
        rectangle = calRectangle_sphere(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
        end
    end

    for atom in atoms
        rectangle = calRectangle_circularThrusters(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end

# 高速化前
function afmize_beta(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.xyz[1, 3*(i-1)+1], tra.xyz[1, 3*(i-1)+2], tra.xyz[1, 3*(i-1)+3],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)

    for h in 1:height, w in 1:width
        probe = Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                      config.range_min.y + (h-0.5) * config.resolution.y,
                      config.probeRadius, config.probeAngle)
        for atom in atoms
            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end

function afmize_gpu(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atom_r = tra.mass
    atom_z = tra.z[:]
    atom_z = atom_z .- (minimum(atom_z) - atom_r[argmin(atom_z)])
    #atom_z .= atom_z .- (minimum(atom_z))
    atom_z_max = maximum(atom_z)
    stage = similar(tra.x, 1, height*width)
    stage .= 0.0

    probe_x = similar(tra.x, 1, width)
    probe_y = similar(tra.x, 1, height)
    probe_x[:] .= range(config.range_min.x, step=config.resolution.x, stop=config.range_min.x + (width-1)*config.resolution.x) .+ (0.5*config.resolution.x)
    probe_y[:] .= range(config.range_min.y, step=config.resolution.y, stop=config.range_min.y + (height-1)*config.resolution.y) .+ (0.5*config.resolution.y)
    probe_r = config.probeRadius
    probe_angle = config.probeAngle

    # collision with sphere
    dist_x2 = (tra.x[:] .- probe_x).^2
    dist_y2 = (tra.y[:] .- probe_y).^2
    index_x_to_xy = [j for j=1:width for i=1:height]
    index_y_to_xy = [i for j=1:width for i=1:height]

    dist_xy2 = dist_x2[:, index_x_to_xy] .+ dist_y2[:, index_y_to_xy]
    #stage_xy = similar(dist_xy2)
    #stage_xy .= 0.0
    dr = probe_r .+ atom_r
    dr2 = dr.^2
    s = dr2 .- dist_xy2
    @show typeof(s)
    index_collide = s .> 0.0
    index_not_collide = .!index_collide
    @show typeof(index_collide)
    if any(index_collide)
        s[Array(index_collide)] .= sqrt.(s[index_collide])
    end
    s .= atom_z .+ s .- probe_r
    if any(index_not_collide)
        s[Array(index_not_collide)] .= 0.0
    end
    s_max = maximum(s, dims=1)
    stage .= max.(stage, s_max)

    # collision with circular thruster
    dist_xy2 .= sqrt.(dist_xy2)
    dist_collision = probe_r .+ atom_r .* cos(probe_angle)
    index_side   = dist_collision .< dist_xy2
    index_corner = (probe_r .< dist_xy2) .& .!index_side
    index_not_collide = .!(index_corner .| index_side)
    s[Array(index_side)] .= (atom_r .* sin.(probe_angle) .- (dist_xy2 .- dist_collision) ./ tan.(probe_angle))[index_side]
    s[Array(index_corner)] .= sqrt.( (atom_r.^2 .- (dist_xy2 .- probe_r).^2)[index_corner]  )
    s .= s .+ atom_z .- probe_r
    s[Array(index_not_collide)] .= 0.0
    s_max = maximum(s, dims=1)
    stage .= max.(stage, s_max)

    stage_reshape = reshape(stage, (height, width))
    return stage_reshape
end

## ---------------------------------------- getafmposteriors_alpha ----------------------------------------------------

function translateafm(afm, (dx, dy))
    afm_translated = zeros(eltype(afm), size(afm))
    (ny, nx) = size(afm)
    for i in maximum([1, 1-dx]):minimum([nx, nx-dx])
        for j in maximum([1, 1-dy]):minimum([ny, ny-dy])
            afm_translated[j+dy, i+dx] = afm[j, i]
        end
    end
    afm_translated
end

function translateafm_periodic(afm, (dx, dy))
    afm_translated = zeros(eltype(afm), size(afm))
    (nx, ny) = size(afm)
    for x in 1:nx
        for y in 1:ny
            xx = x + dx
            yy = y + dy
            if xx > nx
                xx -= nx
            end
            if yy > ny
                yy -= ny
            end
            afm_translated[xx, yy] = afm[x, y]
        end
    end
    afm_translated
end

function direct_convolution(observed, calculated)
    H, W = size(observed)
    ret = zeros(H, W)
    for y in 1:H, x in 1:W
        for dy in 1:H, dx in 1:W
            ny = y + dy - 1
            nx = x + dx - 1
            if ny > H ny -= H end
            if nx > W nx -= W end
            ret[y, x] += observed[dy, dx] * calculated[ny, nx]
        end
    end
    return ret
end

function direct_convolution_non_periodic(observed, calculated)
    H, W = size(observed)
    ret = zeros(H, W)
    for y in 1:H, x in 1:W
        for dy in 1:(H-y+1), dx in 1:(W-x+1)
            ny = y + dy - 1
            nx = x + dx - 1
            ret[y, x] += observed[dy, dx] * calculated[ny, nx]
        end
    end
    return ret
end

function fft_convolution(observed, calculated)
    #return real.(ifft(fft(observed).*conj.(fft(calculated))))
    return real.(ifftshift(ifft(fft(observed).*conj.(fft(calculated)))))
end

function calcLogProb(observed::AbstractMatrix{T}, calculated::AbstractMatrix{T}, convolution_func=fft_convolution) where {T}
    npix = eltype(observed)(size(observed, 1) * size(observed, 2))
    C_o  = sum(observed)
    C_c  = sum(calculated)
    C_cc = sum(calculated.^2)
    C_oo = sum(observed.^2)
    C_oc = convolution_func(observed, calculated)

    log01 = npix .* (C_cc .* C_oo .- C_oc.^2) .+ eltype(observed)(2) .* C_o .* C_oc .* C_c .- C_cc .* C_o.^2 .- C_oo .* C_c.^2
    #log01[Array(log01 .< log(eps(eltype(observed))))] .= eps(eltype(observed))
    log02 = (npix .- eltype(observed)(2)) .* (npix .* C_cc .- C_c.^2)
    #log02 = log02 <= log(eps(eltype(observed))) ? eps(eltype(observed)) : log02
    logprob = eltype(observed)(0.5) .* (eltype(observed)(3.0) .- npix) .* log.(log01) .+ (eltype(observed)(0.5) .* npix .- eltype(observed)(2.0)) .* log(log02)

    return logprob
end

function getafmposterior(afm::Matrix{Float64}, model_array::TrjArray, q_array::Matrix{Float64}, param_array)
    imax_model = 0
    imax_q = 0
    best_model = model_array[1, :]
    best_param = param_array[1]
    best_translate = [0, 0]
    best_afm = similar(afm)
    best_posterior = -Inf

    observed = afm
    npix = size(observed, 1) * size(observed, 2)
    decenter!(model_array)

    posterior_all = zeros(Float64, size(q_array, 1)*length(param_array)*npix, size(model_array, 1))
    posterior = zeros(Float64, size(model_array, 1))

    ### loop over models(structures)
    Threads.@threads for imodel in 1:size(model_array, 1)
        @show imodel
        icount = 0
        model = model_array[imodel, :]
        ### loop over rotations
        for iq in 1:size(q_array, 1)
            q = q_array[iq, :]
            model_rotated = MDToolbox.rotate(model, q)
            ### loop over afmize parameters
            for iparam in 1:length(param_array)
                param = param_array[iparam]
                calculated = afmize(model_rotated, param)
                logprob = calcLogProb(observed, calculated)
                posterior_all[(icount+1):(icount+npix), imodel] .= logprob[:]
                icount += npix
                maximum_logprob = maximum(logprob)
                if best_posterior < maximum_logprob
                    best_posterior = maximum_logprob
                    imax_model = imodel
                    best_model = model_rotated
                    imax_q = iq
                    best_param = param
                    x_center = ceil(Int64, (size(observed,1)/2)+1.0)
                    y_center = ceil(Int64, (size(observed,2)/2)+1.0)
                    dx_estimated = argmax(logprob)[1] - x_center
                    dy_estimated = argmax(logprob)[2] - y_center
                    best_translate = (dx_estimated, dy_estimated)
                    best_afm = translateafm(calculated, best_translate)
                end
            end
        end
    end

    for imodel in 1:size(model_array, 1)
        posterior[imodel] = logsumexp(posterior_all[:, imodel])
    end
    posterior .= posterior .- maximum(posterior)
    posterior .= exp.(posterior)
    posterior .= posterior ./ sum(posterior)

    return imax_model, imax_q, best_param, best_translate, best_afm, best_posterior, posterior
end

function getafmposterior_gpu(afm::AbstractMatrix{T}, model_array::TrjArray{T, U}, q_array::AbstractMatrix{T}, param_array) where {T, U}
    imax_model = 0
    imax_q = 0
    best_param = param_array[1]
    best_translate = [0, 0]
    best_afm = similar(afm)
    best_posterior = -Inf

    observed = afm
    decenter!(model_array)

    for imodel in 1:size(model_array, 1)
        @show imodel
        model = model_array[imodel, :]
        models_rotated = MDToolbox.rotate(model, q_array)
        for iq in 1:size(q_array, 1)
            model_rotated = models_rotated[iq, :]
            for iparam in 1:length(param_array)
                param = param_array[iparam]
                #@show typeof(model_rotated.x)
                calculated = MDToolbox.afmize_gpu(model_rotated, param)
                #@show typeof(calculated)
                logprob = MDToolbox.calcLogProb(observed, calculated, MDToolbox.fft_convolution)
                maximum_logprob = maximum(logprob)
                if best_posterior < maximum_logprob
                    best_posterior = maximum_logprob
                    imax_model = imodel
                    imax_q = iq
                    best_param = param
                    x_center = ceil(Int32, (size(observed,1)/2)+1.0)
                    y_center = ceil(Int32, (size(observed,2)/2)+1.0)
                    dx_estimated = argmax(logprob)[1] - x_center
                    dy_estimated = argmax(logprob)[2] - y_center
                    best_translate = (dx_estimated, dy_estimated)
                    best_afm = translateafm(calculated, best_translate)
                end
            end
        end
    end

    return imax_model, imax_q, best_param, best_translate, best_afm, best_posterior
end

function logprob_eachmodel(model::TrjArray, afm_array, q_array::Matrix{Float64}, param_array)
    nq = size(q_array, 1)
    nparam = length(param_array)
    nafm = length(afm_array)
    npix = size(afm_array[1], 1) * size(afm_array[1], 2)
    logprob_array = zeros(eltype(afm_array[1]), nafm, nq, nparam)
    dxdy_array = Array{Tuple{Int,Int}}(undef, nafm, nq, nparam)
    dxdy = Array{Tuple{Int,Int}}(undef, nafm)
    x_center = ceil(Int64, (size(afm_array[1],2)/2)+1.0)
    y_center = ceil(Int64, (size(afm_array[2],1)/2)+1.0)

    decenter!(model)
    ### loop over rotations
    Threads.@threads for iq in 1:nq
        @show iq
        q = q_array[iq, :]
        model_rotated = MDToolbox.rotate(model, q)
        ### loop over afmize parameters
        for iparam in 1:nparam
            param = param_array[iparam]
            calculated = afmize(model_rotated, param)
            for iafm in 1:nafm
                observed = afm_array[iafm]
                logprob = calcLogProb(observed, calculated)
                logprob_array[iafm, iq, iparam] = logsumexp(logprob)
                dx = argmax(logprob)[2] - x_center
                dy = argmax(logprob)[1] - y_center
                dxdy_array[iafm, iq, iparam] = (dx, dy)
            end
        end
    end

    for iafm in 1:nafm
        ind = argmax(logprob_array[iafm, :, :])
        dxdy[iafm] = dxdy_array[iafm, :, :][ind]
    end

    return (logprob_array=logprob_array, dxdy=dxdy)
end

function getposterior_parallel(models::TrjArray, afm_array, q_array::Matrix{Float64}, param_array)
    nmodel = models.nframe
    nafm = length(afm_array)
    nq = size(q_array, 1)
    nparam = length(param_array)

    p = @showprogress pmap(x -> logprob_eachmodel(x, afm_array, q_array, param_array), models)

    logprob_all = []
    logprob_model = []
    logprob_q = []
    logprob_param = []
    dxdy_best = []
    for iafm = 1:nafm
        pp = zeros(Float64, nmodel, nq, nparam)
        for imodel = 1:nmodel
            for iq = 1:nq
                for iparam = 1:nparam
                    pp[imodel, iq, iparam] = p[imodel].logprob_array[iafm, iq, iparam]
                end
            end
        end

        push!(logprob_all, pp)

        t = zeros(Float64, nmodel)
        for imodel = 1:nmodel
            t[imodel] = logsumexp(pp[imodel, :, :][:])
        end
        push!(logprob_model, t)

        imax = argmax(t)
        push!(dxdy_best, p[imax].dxdy[iafm])

        t = zeros(Float64, nq)
        for iq = 1:nq
            t[iq] = logsumexp(pp[:, iq, :][:])
        end
        push!(logprob_q, t)

        t = zeros(Float64, nparam)
        for iparam = 1:nparam
            t[iparam] = logsumexp(pp[:, :, iparam][:])
        end
        push!(logprob_param, t)
    end

    return (all=logprob_all,
            model=logprob_model,
            q=logprob_q,
            param=logprob_param,
            dxdy=dxdy_best)
end

mutable struct posteriorResult
    each_quate_id
    each_param_id
    posterior_results
    best_translate
    best_posterior
    best_model_id
    best_quate_id
    best_param_id
    best_afm
end

mutable struct posteriorTmpData
    posteriors
    each_best
end

function outputResults(posteriorResults, fileName)
    open(fileName, "w") do out
        N = size(posteriorResults)[1]
        model_size = size(posteriorResults[1].each_quate_id)[1]
        frame_h, frame_w = size(posteriorResults[1].best_afm)
        println(out, "$(N) $(model_size)")
        for result in posteriorResults
            for i in 1:model_size
                if i != 1 print(out, " ") end
                print(out, "$(result.each_quate_id[i])")
            end
            println(out, "")
            for i in 1:model_size
                if i != 1 print(out, " ") end
                print(out, "$(result.each_param_id[i])")
            end
            println(out, "")
            for i in 1:model_size
                if i != 1 print(out, " ") end
                print(out, "$(result.posterior_results[i])")
            end
            println(out, "")
            println(out, "$(result.best_translate[1]) $(result.best_translate[2])")
            println(out, "$(result.best_posterior)")
            println(out, "$(result.best_model_id)")
            println(out, "$(result.best_quate_id)")
            println(out, "$(result.best_param_id)")
            println(out, "$(frame_h) $(frame_w)")
            for h in 1:frame_h
                for w in 1:frame_w
                    if w != 1 print(out, " ") end
                    print(out, "$(result.best_afm[h, w])")
                end
                println(out, "")
            end
        end
    end
end

function inputResults(fileName)
    ret = []
    open(fileName, "r") do io
        N, model_size = map(x->parse(Int,x),split(readline(io)))
        for result in 1:N
            each_quate_id = map(x->parse(Int,x),split(readline(io)))
            each_param_id = map(x->parse(Int,x),split(readline(io)))
            posterior_results = map(x->parse(Float64,x),split(readline(io)))
            best_translate = Tuple(map(x->parse(Int,x),split(readline(io))))
            best_posterior = parse(Float64, readline(io))
            best_model_id = parse(Int, readline(io))
            best_quate_id = parse(Int, readline(io))
            best_param_id = parse(Int, readline(io))
            frame_h, frame_w = map(x->parse(Int,x),split(readline(io)))
            frame = zeros(Float64, frame_h, frame_w)
            for h in 1:frame_h
                frame[h, :] = map(x->parse(Float64,x),split(readline(io)))
            end
            push!(ret, posteriorResult(each_quate_id,
                                        each_param_id,
                                        posterior_results,
                                        best_translate,
                                        best_posterior,
                                        best_model_id,
                                        best_quate_id,
                                        best_param_id,
                                        frame))
        end
    end
    return ret
end

function getafmposteriors_alpha(afm_frames, model_array::TrjArray, quate_array::Matrix{Float64}, param_array, opt = "")
    frame_num = size(afm_frames)[1]
    model_num = size(model_array)[1]
    quate_num = size(quate_array)[1]
    param_num = size(param_array)[1]

    convolution_func = fft_convolution
    if opt == "direct"
        convolution_func = direct_convolution
    end

    println("func is $convolution_func")

    results = []
    tmpData = []
    for i in 1:frame_num
        # 各モデルについて、一つの角度と半径ごとのlogprobの最高値を保持しておく
        posteriors = zeros(model_num, quate_num * param_num)
        # 各モデルについて、最も良い値をだした時のafm画像
        each_quate_id = zeros(Int, model_num)
        each_param_id = zeros(Int, model_num)
        each_best = ones(Float64, model_num) .* -Inf

        best_posterior = -Inf
        best_model_id = 1
        best_quate_id = 1
        best_param_id = 1
        best_afm = zeros(size(afm_frames[1]))
        posterior_results = zeros(model_num)
        push!(results, posteriorResult( each_quate_id,
                                        each_param_id,
                                        posterior_results,
                                        (0, 0),
                                        best_posterior,
                                        best_model_id,
                                        best_quate_id,
                                        best_param_id,
                                        best_afm))
        push!(tmpData, posteriorTmpData(posteriors,
                                        each_best))
    end

    for (model_id, model) in zip(1:model_num, model_array)
        @show model_id
        for quate_id in 1:quate_num
            rotated_model = MDToolbox.rotate(model, quate_array[quate_id, :])
            for param_id in 1:param_num
                cal_frame = afmize(rotated_model, param_array[param_id])
                for frame_id in 1:frame_num
                    prob_mat = calcLogProb(afm_frames[frame_id], cal_frame, convolution_func)
                    max_prob = maximum(prob_mat)
                    id = (quate_id - 1) * param_num + param_id
                    tmpData[frame_id].posteriors[model_id, id] = max_prob

                    if tmpData[frame_id].each_best[model_id] < max_prob
                        tmpData[frame_id].each_best[model_id] = max_prob
                        results[frame_id].each_quate_id[model_id] = quate_id
                        results[frame_id].each_param_id[model_id] = param_id
                    end

                    if results[frame_id].best_posterior < max_prob
                        results[frame_id].best_posterior = max_prob
                        results[frame_id].best_model_id = model_id
                        results[frame_id].best_quate_id = quate_id
                        results[frame_id].best_param_id = param_id
                        results[frame_id].best_translate = Tuple(argmax(prob_mat))
                        results[frame_id].best_afm = translateafm_periodic(cal_frame, results[frame_id].best_translate)
                    end
                end
            end
        end
    end

    for frame_id in 1:frame_num
        tmpData[frame_id].posteriors .-= maximum(tmpData[frame_id].posteriors)
        tmpData[frame_id].posteriors = exp.(tmpData[frame_id].posteriors)
        results[frame_id].posterior_results = sum(tmpData[frame_id].posteriors, dims = 2)
        results[frame_id].posterior_results ./= sum(results[frame_id].posterior_results)
    end

    return results
end


text:
function Base.copy(x::Camera)
    Camera(ntuple(6) do i
        getfield(x, i)
    end...)
end
function Base.:(==)(a::Camera, b::Camera)
    to_value(a.view) == to_value(b.view) &&
    to_value(a.projection) == to_value(b.projection) &&
    to_value(a.resolution) == to_value(b.resolution)
end
function disconnect!(c::Camera)
    for node in c.steering_nodes
        filter!(listeners(node)) do x
            !(x isa CameraLift) # remove all camera lifts
        end
    end
    return
end
function disconnect!(nodes::Vector)
    for node in nodes
        disconnect!(node)
    end
    empty!(nodes)
    return
end

struct CameraLift{F, Args}
    f::F
    args::Args
end

function (cl::CameraLift{F, Args})(val) where {F, Args}
    cl.f(map(to_value, cl.args)...)
end

"""
When mapping over nodes for the camera, we store them in the steering_node vector,
to make it easier to disconnect the camera steering signals later!
"""
function Observables.on(f, c::Camera, nodes::Node...)
    # this basically reimplements onany, which is a bit annoying, but like
    # this we don't have such a closure hell and CameraLift will be nicely
    # identifiable when we disconnect the nodes!
    cl = CameraLift(f, nodes)
    for n in nodes
        on(cl, n)
    end
    push!(c.steering_nodes, nodes...)
    node
end

function Camera(px_area)
    Camera(
        Node(Mat4f0(I)),
        Node(Mat4f0(I)),
        Node(Mat4f0(I)),
        lift(a-> Vec2f0(widths(a)), px_area),
        Node(Vec3f0(1)),
        []
    )
end

function is_mouseinside(scene, target)
    scene === target && return false
    Vec(scene.events.mouseposition[]) in pixelarea(scene)[] || return false
    for child in r.children
        is_mouseinside(child, target) && return true
    end
    return false
end
function is_mouseinside(scene)
    return Vec(scene.events.mouseposition[]) in pixelarea(scene)[]
    # Check that mouse is not inside any other screen
    # for child in scene.children
    #     is_mouseinside(child) && return false
    # end
end

ext: jl


text:
function spharm( l, m, Gvec )

  Ng = size(Gvec)[2]

  const SMALL = 1e-9

  cost = zeros(Float64, Ng)
  sint = zeros(Float64, Ng)
  phi  = zeros(Float64, Ng)

  for ig = 1:Ng
    gmod = sqrt( Gvec[1,ig]^2 + Gvec[2,ig]^2 + Gvec[3,ig]^2 )
    if gmod < SMALL
       cost[ig] = 0.0
    else
       cost[ig] = Gvec[3,ig]/gmod
    end
    #
    #  beware the arc tan, it is defined modulo pi
    #
    if Gvec[1,ig] > SMALL
      phi[ig] = atan( Gvec[2,ig]/Gvec[1,ig] )
    elseif Gvec[1,ig] < -SMALL
      phi[ig] = atan( Gvec[2,ig]/g[1,ig] ) + pi
    else
      phi[ig] = if Gvec[2,ig] >= 0 pi/2 else -pi/2 end
    end
    sint[ig] = sqrt( max( 0.0, 1-cost[ig]^2 ) )
  end

  ylm = zeros( Float64, Ng )

  if l == 0
    for ig = 1:Ng
      ylm[ig] = 0.5*sqrt(1.0/pi)
    end
    return ylm

  elseif l == 1
    # py
    if m == -1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*sint[ig]*sin(phi[ig])
      end
      return ylm
    # pz
    elseif m == 0
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*cost[ig]
      end
      return ylm
    # px
    elseif m == 1
      for ig = 1:Ng
        ylm[ig] = 0.5*sqrt(3/pi)*sint[ig]*cos(phi[ig])
      end
      return ylm
    end

  elseif l == 2
    # dxy
    if m == -2
      for ig = 1:Ng
        #ylm[ig] = sqrt(oppi*15/16)*sin(phi[ig])^2*sin(2*phi[ig])
        ylm[ig] = sqrt(15/16/pi)*sint[ig]^2*sin(2*phi[ig])
      end
      return ylm
    # dyz
    elseif m == -1
      for ig = 1:Ng
        #ylm[ig] = sqrt(oppi*15/4)*cost[ig]*sint[ig]*sin(phi[ig])
        ylm[ig] = sqrt(15/(4*pi))*cost[ig]*sint[ig]*sin(phi[ig])
      end
      return ylm
    # dz2
    elseif m == 0
      for ig = 1:Ng
        #ylm[ig] = sqrt(oppi*5/8)*( 2*cost[ig]^2 - sin(phi[ig])^2 )
        ylm[ig] = 0.25*sqrt(5/pi)*( 3*cost[ig]^2 - 1 )
      end
      return ylm
    # dxz
    elseif m == 1
      for ig = 1:Ng
        #ylm[ig] = sqrt(oppi*15/4)*cost[ig]*sint[ig]*cos(phi[ig])
        ylm[ig] = sqrt(15/4/pi)*cost[ig]*sint[ig]*cos(phi[ig])
      end
      return ylm
    # dx2-y2
    elseif m == 2
      for ig = 1:Ng
        #ylm[ig] = sqrt(oppi*15/16)*sin(phi[ig])^2*cos(2*phi[ig])
        ylm[ig] = 0.5*sqrt(15/4/pi)*sint[ig]^2*cos(2*phi[ig])
      end
      return ylm

    end

  end

end


text:
# see good tutorial here: http://ksimek.github.io/2013/08/13/intrinsic/
# the camera here is initialized to look towards -z, and up is y, just like the default Blender camera

function _getIntrinsic(FOVH)
    # focal length, asumming some detector size which is irrelevant
    fl = 1.0/tan(FOVH/2.0)
    [fl   0     0   0;
      0  fl     0   0;
      0   0     1   0;
      0   0     0   1]
end

function _getExtrinsic(ROT_XYZ, TRANS_XYZ)
    R = (rotXYZ(ROT_XYZ) * rotXYZ([pi,0.0,pi]))'
    ex = zeros(4,4)
    ex[1:3,1:3] = R
    ex[1:3,4] = -R * TRANS_XYZ
    ex[4,4] = 1.0
    ex
end

function rotX(theta)
    ct = cos(theta)
    st = sin(theta)
    [1   0    0;
     0   ct  -st;
     0   st   ct]
end

function rotY(theta)
    ct = cos(theta)
    st = sin(theta)
    [ ct   0    st;
      0    1     0;
     -st   0    ct]
end

function rotZ(theta)
    ct = cos(theta)
    st = sin(theta)
    [ct   -st    0;
     st    ct    0;
     0     0     1]
end

rotXYZ(XYZ::Vector{Float64}) = rotZ(XYZ[3]) * rotY(XYZ[2]) * rotX(XYZ[1])

Camera(ROT_XYZ::Vector{Float64}, TRANS_XYZ::Vector{Float64}, FOVH::Float64) =
    _getIntrinsic(FOVH) * _getExtrinsic(ROT_XYZ, TRANS_XYZ)

@inline function xyz_to_image_coords(CAM::Array{Float64,2}, P_XYZ::Vector{Float64}, image_width::Int64, image_height::Int64)
    i = CAM * [P_XYZ..., 1.0]
    imx = i[1] / i[3] * image_width / 2+image_width/2
    imy = i[2] / i[3] * image_height / 2*image_width/image_height+image_height/2
    image_width-imx, image_height-imy, i[3]
end

function line_to_2d(CAM, xl,yl,W=256,H=192)
    x2d_l = zeros(length(xl))
    y2d_l = zeros(length(yl))
    z2d_l = zeros(length(yl))
    for i in eachindex(xl)
        if isnan(xl[i])
            x2d_l[i] = NaN
            y2d_l[i] = NaN
            z2d_l[i] = NaN
            continue
        end
        v = xyz_to_image_coords(CAM, [xl[i],yl[i],0.0], W, H)
        x2d_l[i] = v[1]
        y2d_l[i] = v[2]
        z2d_l[i] = v[3]
    end
    ixs = find((z2d_l.>0.0) .| (isnan.(z2d_l)))
    x2d_l[ixs],y2d_l[ixs]
end


text:
@testset "file2 gets run" begin
    global file2_run
    @test !file2_run
    file2_run = true
    @test file2_run
end


text:
module GS

include("fea.jl")
include("otm.jl")
include("graficos.jl")

using .FEA,
      .Otm,
      .Graficos

# FEA
export No,
       Material,
       Elemento,
       Estrutura,
       comprimento,
       num_gls_estrutura,
       vetor_forcas,
       volume,
       deslocamentos,
       gls,
       k,
       tem_sobreposicao,
       carregar_apoiar_no!

# Otm
export Otimizador,
       otimizar!,
       gerar_trelica_hiperconectada

# Graficos
export plotar_estrutura_inicial,
       plotar_estrutura_otimizada

end

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
#
# Iterative Halo orbit and manifold solvers
#

"""
Returns the Monodromy Matrix for a Halo orbit.
"""
function monodromy(orbit::CircularRestrictedThreeBodyOrbit, T; reltol = 1e-14, abstol = 1e-14)
    orb = NormalizedSynodicSTMCR3BPOrbit(normalize(synodic(orbit)))
    problem = ODEProblem(orb, T)

    final = solve(problem; reltol=reltol, abstol=abstol, save_everystep=false).u[end]

    SMatrix{6,6}(final[7:end]...) |> Matrix
end

"""
Returns true if a `RestrictedThreeBodySystem` is numerically periodic.
"""
function isperiodic(orbit::CircularRestrictedThreeBodyOrbit, T; reltol = 1e-14, abstol = 1e-14, atol = 1e-6)
    orb = NormalizedSynodicSTMCR3BPOrbit(normalize(synodic(orbit)))
    problem = ODEProblem(orbit, T)

    final = solve(problem; reltol=reltol, abstol=abstol, save_everystep=false).u[end]

    return all(isapprox.(position_vector(orbit.state), final[1:3]; atol = atol)) && all(isapprox.(velocity_vector(orbit.state), final[4:6]; atol = atol))
end

"""
Returns a numerical solution for a Halo orbit about `L`.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `Az`: Desired non-dimensional Z-amplitude for Halo orbit.
- `ϕ`: Desired Halo orbit phase.
- `L`: Lagrange point to orbit (L1 or L2).
- `hemisphere`: Specifies northern or southern Halo orbit.

__Outputs:__
- Tuple of initial states: `(r, v)` where `r::Vector{<:AbstractFloat}`, `v::Vector{<:Abstractfloat}`.
- Throws `ArgumentError` if L is not `:L1` or `:L2`

__References:__

The iterative scheme was pulled from directly from literature
and sample code, including Rund 2018,
and Dr. Mireles' lecture notes and EarthSunHaloOrbit_NewtonMewhod.m
file available on their website.
Specifically, the half-period iterative scheme (the `F` matrix
in the source code, and corresponding "next guess" computation)
was ported __directly__ from Dr. Mireles' public code and notes, which are
available online.
- [Dr. Mireles Notes](http://cosweb1.fau.edu/~jmirelesjames/hw5Notes.pdf)
- [Dr. Mireles Code](http://cosweb1.fau.edu/~jmirelesjames/notes.html)
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/).
"""
function halo(μ; Az=0.0, L=1, hemisphere=:northern,
              tolerance=1e-8, max_iter=100,
              reltol=1e-14, abstol=1e-14,
              nan_on_fail = true, disable_warnings = false)

    r₀, v₀, Τ = analyticalhalo(μ; Az=Az, ϕ=0.0, L=L, hemisphere=hemisphere)
    r₀ = r₀[1,:]
    v₀ = v₀[1,:]
    τ  = Τ/2
    δu = Vector{typeof(τ)}(undef, 6)
    Id = Matrix(I(6))

    for i ∈ 1:max_iter

        problem = ODEProblem(
            AstrodynamicalModels.CR3BPWithSTMVectorField,
            vcat(r₀, v₀, Id...),
            (0.0, τ),
            (μ,)
        )

        retcode, rₛ, vₛ, Φ = let
            sols  = solve(problem, Vern9(); reltol=reltol, abstol=abstol)
            final = sols.u[end]
            sols.retcode, final[1:3], final[4:6], SMatrix{6,6}(final[7:end]...)
        end

        AstrodynamicalModels.CR3BPVectorField(δu, vcat(rₛ, vₛ), (μ,), NaN)
        ∂vₛ = δu[4:6]

        # All code in this `if, else, end` block is ported from
        # Dr. Mireles' MATLAB code, which is available on his
        # website: http://cosweb1.fau.edu/~jmirelesjames/notes.html.
        # Lecture notes, which describe this algorithm further,
        # are available for reference as well:
        # http://cosweb1.fau.edu/~jmirelesjames/hw5Notes.pdf
        if Az ≉ 0
            F = @SMatrix [
                Φ[4,1] Φ[4,5] ∂vₛ[1];
                Φ[6,1] Φ[6,5] ∂vₛ[3];
                Φ[2,1] Φ[2,5]  vₛ[2]
            ]

            xᵪ = SVector(r₀[1], v₀[2], τ) - inv(F) * SVector(vₛ[1], vₛ[3], rₛ[2])

            r₀[1] = xᵪ[1]
            v₀[2] = xᵪ[2]
            τ     = xᵪ[3]
        else
            F = @SMatrix [
                Φ[4,3] Φ[4,5] ∂vₛ[1];
                Φ[6,3] Φ[6,5] ∂vₛ[3];
                Φ[2,3] Φ[2,5]  vₛ[2]
            ]

            xᵪ = SVector(r₀[3], v₀[2], τ) - inv(F) * SVector(vₛ[1], vₛ[3], rₛ[2])

            r₀[3] = xᵪ[1]
            v₀[2] = xᵪ[2]
            τ     = xᵪ[3]
        end

        if abs(vₛ[1]) ≤ tolerance && abs(vₛ[3]) ≤ tolerance
            break;
        elseif τ > 5 * one(τ)
            disable_warnings || @warn "Unreasonably large halo period, $τ, ending iterations."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        elseif i == max_iter
            disable_warnings || @warn "Desired tolerance was not reached, and iterations have hit the maximum number of iterations: $max_iter."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        elseif retcode != :Success
            !disable_warnings || @warn "Integrator returned $(string(retcode))."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        end
    end

    return r₀, v₀, 2τ

end

"""
A `halo` wrapper! Returns a `CircularRestrictedThreeBodyOrbit`.
Returns a tuple: `halo_orbit, halo_period`.
"""
function halo(sys::CircularRestrictedThreeBodySystem;
              Az=0.0,
              L=1,
              hemisphere=:northern,
              tolerance=1e-8,
              max_iter=20,
              reltol=1e-14,
              abstol=1e-14,
              nan_on_fail = true,
              disable_warnings = false)

    if Az isa Unitful.Length
        Az = nondimensionalize_length(Az, normalized_length_unit(sys))
    end

    r,v,T = halo(normalized_mass_parameter(sys);
                 Az = Az, L = L, hemisphere = hemisphere,
                 tolerance = tolerance, max_iter = max_iter, reltol = reltol, abstol = abstol,
                 nan_on_fail = nan_on_fail, disable_warnings = disable_warnings)
    orbit = CircularRestrictedThreeBodyOrbit(redimensionalize_length.(r, normalized_length_unit(sys)),
                                             redimensionalize_velocity.(v, normalized_length_unit(sys), normalized_time_unit(sys)),
                                             sys) |> normalize
    return orbit, T
end

"""
Calculates the eigenvector associated with the __stable manifold__
of a Monodromy matrix.
"""
function stable_eigenvector(monodromy::AbstractMatrix; check=true)
    evals, evecs = eigen(monodromy)
    evals = filter(isreal, evals) .|> real
    evecs = filter(x->!isempty(x), map(vec->filter(x->all(isreal.(vec)), vec), eachcol(evecs))) .|> real

    imin = findmin(evals)[2]
    imax = findmax(evals)[2]

    !check || @assert (evals[imin] * evals[imax]) ≈ one(eltype(evals)) "Min and max eigenvalue should be multiplicative inverses. Invalid Monodromy matrix. Product equals $(evals[imin] * evals[imax]), not $(one(eltype(evals)))."
    return evecs[imin]
end

"""
Calculates the eigenvector associated with the __unstable manifold__
of a Monodromy matrix.
"""
function unstable_eigenvector(monodromy::AbstractMatrix; check=true)
    evals, evecs = eigen(monodromy)
    evals = filter(isreal, evals) .|> real
    evecs = filter(x->!isempty(x), map(vec->filter(x->all(isreal.(vec)), vec), eachcol(evecs))) .|> real

    imin = findmin(evals)[2]
    imax = findmax(evals)[2]

    !check || @assert (evals[imin] * evals[imax]) ≈ one(eltype(evals)) "Min and max eigenvalue should be multiplicative inverses. Invalid Monodromy matrix. Product equals $(evals[imin] * evals[imax]), not $(one(eltype(evals)))."
    return evecs[imax]
end

"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a nearby `CircularRestrictedThreeBodyOrbit` on the
__stable manifold__ of the initial state.
"""
function manifold(orbit::NormalizedSynodicSTMCR3BPOrbit, V::AbstractVector; eps = 1e-8)
    V = normalize(V)
    r = orbit.state.cart.r + eps * V[1:3]
    v = orbit.state.cart.v + eps * V[4:6]

    cart = CartesianState(r, v, orbit.state.cart.t, Synodic;
                          lengthunit = NormalizedLengthUnit(),
                          timeunit = NormalizedTimeUnit())

    return CircularRestrictedThreeBodyOrbit(cart, orbit.system)
end

"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a collection of `Trajectory` instances representing
the stable manifold near orbit.
"""
function stable_manifold(orbit::NormalizedSynodicCR3BPOrbit, T::Real;
                         duration = T, reltol = 1e-14,
                         abstol = 1e-14, eps = 1e-8,
                         saveat = 1e-2, num_trajectories = 10, parallel = false)

    @assert duration > 0 "The provided duration cannot be zero or negative. Duration provided was $duration"
    @assert isperiodic(orbit, T) "The provided orbit is not periodic!"

    M = monodromy(orbit, T)
    V = stable_eigenvector(M)

    orbit = NormalizedSynodicSTMCR3BPOrbit(orbit)
    traj = propagate(orbit, T; reltol = reltol, abstol = abstol, saveat = saveat)
    if length(traj) ÷ num_trajectories < 1
        ind = 1:length(traj)
    else
        ind  = 1 : length(traj) ÷ num_trajectories : length(traj) - (length(traj) % num_trajectories)
    end

    f = step -> reverse(propagate(manifold(step, step.state.stm * V; eps = eps), -duration;
                                  reltol = reltol, abstol = abstol))
    return parallel ? pmap(f, traj[ind]) : map(f, traj[ind])
end

"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a collection of `Trajectory` instances representing
the unstable manifold near orbit.
"""
function unstable_manifold(orbit::NormalizedSynodicCR3BPOrbit, T::Real;
                           duration = T, reltol = 1e-14,
                           abstol = 1e-14, eps = 1e-8,
                           num_trajectories = 10,
                           saveat = 1e-2, parallel = false)

    @assert duration > zero(duration) "The provided duration cannot be zero or negative."
    @assert isperiodic(orbit, T) "The provided orbit is not periodic!"

    orbit = NormalizedSynodicSTMCR3BPOrbit(orbit)
    M = monodromy(orbit, T)
    V = unstable_eigenvector(M)

    traj = propagate(orbit, T; reltol = reltol, abstol = abstol, saveat = saveat)
    if length(traj) ÷ num_trajectories < 1
        ind = 1:length(traj)
    else
        ind  = 1 : length(traj) ÷ num_trajectories : length(traj) - (length(traj) % num_trajectories)
    end

    f = step -> propagate(manifold(step, step.state.stm * V; eps = eps), duration;
                          reltol = reltol, abstol = abstol, saveat = saveat)
    return parallel ? pmap(f, traj[ind]) : map(f, traj[ind])
end


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

using Main.aimajulia.utils;

using LinearAlgebra;

#The following learning tests are from the aima-python doctests

@test (repr(euclidean_distance([1, 2], [3, 4])) == "2.8284271247461903");

@test (repr(euclidean_distance([1, 2, 3], [4, 5, 6])) == "5.196152422706632");

@test (repr(euclidean_distance([0, 0, 0], [0, 0, 0])) == "0.0");

@test (root_mean_square_error([2, 2], [2, 2]) == 0);

@test (root_mean_square_error([0, 0], [0, 1]) == sqrt(0.5));

@test (root_mean_square_error([1, 0], [0, 1]) == 1);

@test (root_mean_square_error([0, 0], [0, -1]) == sqrt(0.5));

@test (root_mean_square_error([0, 0.5], [0, -0.5]) == sqrt(0.5));

@test (manhattan_distance([2, 2], [2, 2]) == 0);

@test (manhattan_distance([0, 0], [0, 1]) == 1);

@test (manhattan_distance([1, 0], [0, 1]) == 2);

@test (manhattan_distance([0, 0], [0, -1]) == 1);

@test (manhattan_distance([0, 0.5], [0, -0.5]) == 1);

@test (mean_boolean_error([1, 1], [0, 0]) == 1)

@test (mean_boolean_error([0, 1], [1, 0]) == 1)

@test (mean_boolean_error([1, 1], [0, 1]) == 0.5)

@test (mean_boolean_error([0, 0], [0, 0]) == 0)

@test (mean_boolean_error([1, 1], [1, 1]) == 0)

@test (mean_error([2, 2], [2, 2]) == 0);

@test (mean_error([0, 0], [0, 1]) == 0.5);

@test (mean_error([1, 0], [0, 1]) ==  1);

@test (mean_error([0, 0], [0, -1]) ==  0.5);

@test (mean_error([0, 0.5], [0, -0.5]) == 0.5);

@test (gaussian(1,0.5,0.7) == 0.6664492057835993);

@test (gaussian(5,2,4.5) == 0.19333405840142462);

@test (gaussian(3,1,3) == 0.3989422804014327);

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv", exclude=[4]);

@test (iris_dataset.inputs == [1, 2, 3]);

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");
means_dict, deviations_dict = find_means_and_deviations(iris_dataset);

@test (means_dict["setosa"][1] == 5.006);

@test (means_dict["versicolor"][1] == 5.935999999999999);

@test (means_dict["virginica"][1] == 6.587999999999999);

@test (deviations_dict["setosa"][1] == 0.3524896872134513);

@test (deviations_dict["versicolor"][1] == 0.5161711470638634);

@test (deviations_dict["virginica"][1] == 0.6358795932744321);

cpd = CountingProbabilityDistribution();

for i in 1:10000
    add(cpd, rand(RandomDeviceInstance, ["1", "2", "3", "4", "5", "6"]));
end

probabilities = collect(cpd[n] for n in ("1", "2", "3", "4", "5", "6"));

@test ((1.0/7.0) <= reduce(min, probabilities) <= reduce(max, probabilities) <= (1.0/5.0));

zoo_dataset = DataSet(name="zoo", examples="./aima-data/zoo.csv");

pl = PluralityLearner(zoo_dataset);

@test (predict(pl, [1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 4, 1, 0, 1]) == "mammal");

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");

# Naive Discrete Model
nbdm = NaiveBayesLearner(iris_dataset, continuous=false);

@test (predict(nbdm, [5, 3, 1, 0.1]) == "setosa");

@test (predict(nbdm, [6, 3, 4, 1.1]) == "versicolor");

@test (predict(nbdm, [7.7, 3, 6, 2]) == "virginica");

# Naive Continuous Model
nbcm = NaiveBayesLearner(iris_dataset, continuous=true);

@test (predict(nbcm, [5, 3, 1, 0.1]) == "setosa");

@test (predict(nbcm, [6, 5, 3, 1.5]) == "versicolor");

@test (predict(nbcm, [7, 3, 6.5, 2]) == "virginica");

# Naive Conditional Probability Model
d1 = CountingProbabilityDistribution(vcat(fill('a', 50), fill('b', 30), fill('c', 15)));
d2 = CountingProbabilityDistribution(vcat(fill('a', 30), fill('b', 45), fill('c', 20)));
d3 = CountingProbabilityDistribution(vcat(fill('a', 20), fill('b', 20), fill('c', 35)));
nbsm = NaiveBayesLearner(Dict([(("First", 0.5), d1), (("Second", 0.3), d2), (("Third", 0.2), d3)]), simple=true);

@test (predict(nbsm, "aab") == "First");

@test (predict(nbsm, ['b', 'b']) == "Second");

@test (predict(nbsm, "ccbcc") == "Third");

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");

k_nearest_neighbors = NearestNeighborLearner(iris_dataset, 3);

@test (predict(k_nearest_neighbors, [5, 3, 1, 0.1]) == "setosa");

@test (predict(k_nearest_neighbors, [6, 5, 3, 1.5]) == "versicolor");

@test (predict(k_nearest_neighbors, [7.5, 4, 6, 2]) == "virginica");

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");

# Truncated SVD
test_epsilon = 1e-9

# Test for consistent properties of the SVD, recall that we randomize initial values during computation
for test_iteration in 1:1000
    U, S, V = svd([1 5 9 13; 2 6 10 14; 3 7 11 15; 4 8 12 16]);
    ev_m, ev_n, singular_values = truncated_svd([1 5 9 13; 2 6 10 14; 3 7 11 15; 4 8 12 16]);

    @test (all(i->(abs(i) < test_epsilon), singular_values - S[1:2]));

    # A * v = eigenvalue * v
    @test (all(i->(abs(i) < test_epsilon),
                (ev_n[1]*singular_values[1])
                - [1 5 9 13; 2 6 10 14; 3 7 11 15; 4 8 12 16]*ev_m[1]));
    @test (all(i->(abs(i) < test_epsilon),
                (ev_n[2]*singular_values[2])
                - [1 5 9 13; 2 6 10 14; 3 7 11 15; 4 8 12 16]*ev_m[2]));

    # check if returned eigenvectors applies to the transpose of our original matrix
    @test (all(i->(abs(i) < test_epsilon),
                (ev_m[1]*singular_values[1])
                - [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]*ev_n[1]));
    @test (all(i->(abs(i) < test_epsilon),
                (ev_m[2]*singular_values[2])
                - [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]*ev_n[2]));

    #A=U*S*V^T
    @test (all(i->(abs(i) < test_epsilon),
                ([1.0 5.0 9.0 13.0; 2.0 6.0 10.0 14.0; 3.0 7.0 11.0 15.0; 4.0 8.0 12.0 16.0]
                - (hcat(ev_n[1], ev_n[2]) * Diagonal(singular_values) * transpose(hcat(ev_m[1], ev_m[2]))))));
end

# aima-python doctests for truncated_svd()
@test (all(i->(abs(i) < test_epsilon),
            (truncated_svd([17 0; 0 11])[3] - [17.0, 11.0])));

@test (all(i->(abs(i) < test_epsilon),
            (truncated_svd([17 0; 0 -34])[3] - [34.0, 17.0])));

@test (all(i->(abs(i) < test_epsilon),
            (truncated_svd([1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0])[3] - [3.0, sqrt(5)])));

@test (all(i->(abs(i) < test_epsilon),
            (truncated_svd([3 2 2; 2 3 -2])[3] - [5.0, 3.0])));

# Decision Tree
dtl = DecisionTreeLearner(iris_dataset);

@test (predict(dtl, [5, 3, 1, 0.1]) == "setosa");

@test (predict(dtl, [6, 5, 3, 1.5]) == "versicolor");

@test (predict(dtl, [7.5, 4, 6, 2]) == "virginica");

function test_rf_predictions(ex1_results::AbstractVector, ex2_results::AbstractVector, ex3_results::AbstractVector)
    local rf::RandomForest = RandomForest(iris_dataset);
    push!(ex1_results, (predict(rf, [5, 3, 1, 0.1]) == "setosa"));
    push!(ex2_results, (predict(rf, [6, 5, 3, 1]) == "versicolor"));
    push!(ex3_results, (predict(rf, [7.5, 4, 6, 2]) == "virginica"));
    nothing;
end

setosa_results = Array{Bool, 1}();
versicolor_results = Array{Bool, 1}();
virginica_results = Array{Bool, 1}();
# Run test_rf_predictions() 1000 times.
println("@time for i in 1:1000\n\ttest_rf_predictions(setosa_results, versicolor_results, virginica_results);\nend");
@time for i in 1:1000
    test_rf_predictions(setosa_results, versicolor_results, virginica_results);
end

setosa_results_count = count((function(b::Bool)
                                    return b;
                                end), setosa_results);
versicolor_results_count = count((function(b::Bool)
                                    return b;
                                end), versicolor_results);
virginica_results_count = count((function(b::Bool)
                                    return b;
                                end), virginica_results);

# lowest setosa_results_count result previously obtained was 970
@test (setosa_results_count >= 960);

# lowest versicolor_results_count result previously obtained was 959
@test (versicolor_results_count >= 950);

# lowest virginica_results_count result previously obtained was 996
@test (virginica_results_count >= 990);

println();
println("setosa assert count (out of 1000): ", setosa_results_count);
println("setosa assertion failure rate: approximately ", Float64(1000 - setosa_results_count)/10.0, "%");
println("versicolor assert count (out of 1000): ", versicolor_results_count);
println("versicolor assertion failure rate: approximately ", Float64(1000 - versicolor_results_count)/10.0, "%");
println("virginica assert count (out of 1000): ", virginica_results_count);
println("virginica assertion failure rate: approximately ", Float64(1000 - virginica_results_count)/10.0, "%");
println();

weights = random_weights(-0.5, 0.5, 10);

@test (length(weights) == 10);

@test (all(((weight >= -0.5) && (weight <= 0.5)) for weight in weights));

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");

# The DataType of the example classification must match the eltype of the classes array.

classes = map(SubString{String}, ["setosa", "versicolor", "virginica"]);

classes_to_numbers(iris_dataset, classes);

nnl = NeuralNetworkLearner(iris_dataset, hidden_layers_sizes=[5], learning_rate=0.15, epochs=75);

neural_network_learner_score = aimajulia.grade_learner(nnl,
                                                        [([5, 3, 1, 0.1], 1),
                                                        ([5, 3.5, 1, 0], 1),
                                                        ([6, 3, 4, 1.1], 2),
                                                        ([6, 2, 3.5, 1], 2),
                                                        ([7.5, 4, 6, 2], 3),
                                                        ([7, 3, 6, 2.5], 3)]);

println("neural network learner score (out of 1.0): ", neural_network_learner_score);

# Allow up to 2 failed tests.
@test (neural_network_learner_score >= (2/3));

neural_network_learner_error_ratio = aimajulia.error_ratio(nnl, iris_dataset);

println("neural network learner error ratio: ", (neural_network_learner_error_ratio * 100), "%");
println();

# NeuralNetworkLearner previously had an error ratio of 0.33333333333333337.
@test (neural_network_learner_error_ratio < 0.40);

iris_dataset = DataSet(name="iris", examples="./aima-data/iris.csv");

classes_to_numbers(iris_dataset, nothing);

pl = PerceptronLearner(iris_dataset);

perceptron_learner_score = aimajulia.grade_learner(pl,
                                                    [([5, 3, 1, 0.1], 1),
                                                    ([5, 3.5, 1, 0], 1),
                                                    ([6, 3, 4, 1.1], 2),
                                                    ([6, 2, 3.5, 1], 2),
                                                    ([7.5, 4, 6, 2], 3),
                                                    ([7, 3, 6, 2.5], 3)]);

println("perceptron learner score (out of 1.0): ", perceptron_learner_score);

# Allow up to 3 failed tests.
@test (perceptron_learner_score >= (1/2));

perceptron_learner_error_ratio = aimajulia.error_ratio(pl, iris_dataset);

println("perceptron learner error ratio: ", (perceptron_learner_error_ratio * 100), "%");
println();

@test (perceptron_learner_error_ratio < 0.40);

@test (weighted_mode("abbaa", [1, 2, 3, 1, 2]) == "b");

@test (weighted_mode(["a", "b", "b", "a", "a"], [1, 2, 3, 1, 2]) == "b");

@test (weighted_replicate(["A", "B", "C"], [1, 2, 1], 4) == ["A", "B", "B", "C"]);

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork