text:
using Distributed
@everywhere using Random
@everywhere using Statistics
@everywhere using LinearAlgebra  # for using eigen
using Plots
using StatsPlots
using SharedArrays
using StatsBase     # for fitting histograms
using CSV, DataFrames   # for saving results as CSV for later
using LaTeXStrings

##### Some local setups
# If run from command line externally:
#projDir = dirname(@__FILE__)

# If run from Jupyter/Hydrogen, maybe change to suit you:
projDir= "/lhome/johhub/Desktop/ABDA/A7"

#####
#%% Input Data #################################################################
@everywhere y = [607, 583, 521, 494, 369, 782, 570, 678, 467, 620, 425, 395, 346, 361, 310,
300, 382, 294, 315, 323, 421, 339, 398, 328, 335, 291, 329, 310, 294, 321, 286,
349, 279, 268, 293, 310, 259, 241, 243, 272, 247, 275, 220, 245, 268, 357, 273,
301, 322, 276, 401, 368, 149, 507, 411, 362, 358, 355, 362, 324, 332, 268, 259,
274, 248, 254, 242, 286, 276, 237, 259, 251, 239, 247, 260, 237, 206, 242, 361,
267, 245, 331, 357, 284, 263, 244, 317, 225, 254, 253, 251, 314, 239, 248, 250,
200, 256, 233, 427, 391, 331, 395, 337, 392, 352, 381, 330, 368, 381, 316, 335,
316, 302, 375, 361, 330, 351, 186, 221, 278, 244, 218, 126, 269, 238, 194, 384,
154, 555, 387, 317, 365, 357, 390, 320, 316, 297, 354, 266, 279, 327, 285, 258,
267, 226, 237, 264, 510, 490, 458, 425, 522, 927, 555, 550, 516, 548, 560, 545,
633, 496, 498, 223, 222, 309, 244, 207, 258, 255, 281, 258, 226, 257, 263, 266,
238, 249, 340, 247, 216, 241, 239, 226, 273, 235, 251, 290, 473, 416, 451, 475,
406, 349, 401, 334, 446, 401, 252, 266, 210, 228, 250, 265, 236, 289, 244, 327,
274, 223, 327, 307, 338, 345, 381, 369, 445, 296, 303, 326, 321, 309, 307, 319,
288, 299, 284, 278, 310, 282, 275, 372, 295, 306, 303, 285, 316, 294, 284, 324,
264, 278, 369, 254, 306, 237, 439, 287, 285, 261, 299, 311, 265, 292, 282, 271,
268, 270, 259, 269, 249, 261, 425, 291, 291, 441, 222, 347, 244, 232, 272, 264,
190, 219, 317, 232, 256, 185, 210, 213, 202, 226, 250, 238, 252, 233, 221, 220,
287, 267, 264, 273, 304, 294, 236, 200, 219, 276, 287, 365, 438, 420, 396, 359,
405, 397, 383, 360, 387, 429, 358, 459, 371, 368, 452, 358, 371]

@everywhere ind = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8,
 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
  11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
   12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14,
   14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16,
   16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20,
   20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21,
   21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,
   22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24,
   24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25,
   25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28,
   28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30,
   30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34,
   34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34];

# Indicator for each individual j whether he/she is a child or not:
@everywhere child_j = [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1,
           0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0];

# Indicator for each observation i whether it comes from a child or not:
@everywhere child_i = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];

# Attempt numbers of each observation i:
@everywhere x = [1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3,
4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 17, 18, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4,
5, 6, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, 7,
8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 1, 2, 3, 4,
5, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18];

# Alternatively, have it all in one array:
@everywhere x_alt = Array{Int64,1}[
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6, 7, 8, 9],
[1], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13],
[1, 2, 3, 4, 5, 6],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6, 7, 8, 9],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[1, 2, 3, 4, 5, 6],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]];

#####
#%% Functions ##################################################################
# @everywhere to make them known to all workers
# Progress bar from Jesper:
@everywhere function progressbar(i,N, msg = "")
    if mod(i,round(N/1000))==0
        print("\e[2K") # clear whole line
        print("\e[1G") # move cursor to column 1
        print(msg * " => " * string(round(i/N*100,sigdigits=4)) * "%")
    end
end

@everywhere function mySliceSampler(pdf_log, x_0, w, N = 2000; burnIn = 1000, m=100, verbose=false, msg="")
    # REQUIRES MODULES:
    # Random
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_0 --> inital vector of x (dimension Dx1 or 1xD), where D dimension of pdf
    # w --> typical slice size, for a window (dimension Dx1 or 1xD)
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    # burnIn --> number of burn-in samples (optional, results will be discarded)
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values (Nx1)

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations
    pdflog_x_1 = 0.0;

    # to be able to input scalars:
    #if typeof(x_0) == Float64
    #    x_0 = [x_0]
    #end

    for ii in 1:(N+burnIn)
        if verbose
            progressbar(ii,N+burnIn,msg)
        end

        # Update the new x progressively for each dimension:
        L   = 1*x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   = 1*x_0
        x_1 = 1*x_0  # "1*" needed, for correct type
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            if (ii == 1)
                #vertical = rand() * pdf_log(x_0) # in normal space
                z = pdf_log(x_0) + log(rand()) # in log-space
            else
                # saving evaluations, pdflog_x_1 was updated for current permutation in the end
                z = pdflog_x_1 + log(rand())
            end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            L[d] = x_0[d] - rand() * w[d]
            R[d] = L[d] + w[d]

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L[d] -= w[d]
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R[d] += w[d]
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                x_1[d] = L[d] + rand() * (R[d] - L[d])

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0[d] = x_1[d]     # update value for this dimension
                    break
                # Value was not within slice, shrink the interval:
                elseif (x_1[d] < x_0[d])
                    L[d] = x_1[d]
                elseif (x_1[d] > x_0[d])
                    R[d] = x_1[d]
                else
                    throw(ErrorException("Error during shrinking the interval"))
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        # Just overwrite the burnIn by using modular arithmetic:
        i = 1 + (ii-1)%N    # gives 1,2,3,...N,1,2,3,...,(N+burnIn)
        x_s[i,:] = x_1     # different directions in arrays n'importe, sample point, after one round of permutations
        pdflog_x_s[i] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations
    end
    return x_s, pdflog_x_s
end

@everywhere function mySliceSampler_Cov(pdf_log, x_0, C, N; burnIn = 1000, m=100, verbose=true, msg="")
    # REQUIRES MODULES:
    # Random
    # LinearAlgebra
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_0 --> inital vector (mode of a previous run) of x (dimension 1xD), where D dimension of pdf
    # C --> Covariance matrix of a previous burn-in of the same log-pdf
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    # burnIn --> number of burn-in samples (optional, results will be discarded)
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values

    λ,Eigvec = LinearAlgebra.eigen(C)        # eigenvalue and matrix of eigenvectors, C is the covariance matrix of the coordinates of a previous run
    w = sqrt.(abs.(λ))    # the single standard-deviations along the dimensions

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations

    pdflog_x_1 = pdf_log(x_0);       # first value is already known!
    wd = zeros(D)       # preallocating window

    L = similar(x_0)
    R = similar(x_0)
    x_1 = similar(x_0)
    # Using dots for array assignment is a bit faster

    for ii in 1:(N+burnIn)
        if verbose
            progressbar(ii,N+burnIn,msg)
        end

        # Update the new x progressively for each dimension:
        L   .= x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   .= x_0
        x_1 .= x_0  # "1*" needed, for correct type
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            #vertical = rand() * pdf_log(x_0) # in normal space
            z = pdflog_x_1 + log(rand())                # saving evaluations, pdflog_x_1 was updated for current permutation in the end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            # Here the entire boundaries are updates simultaneously, to
            # account for the covariation between parameters!
            wd = w[d]*Eigvec[:,d]    # associated co-std value times its eigenvector gives estimate of window
            L .= x_0 .- rand() .* wd
            R .= L .+ wd

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L .-= wd
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R .+= wd
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                # Updating over all dimensions simultaneously, to account for the
                # covariation between parameters:
                x_1 .= L .+ rand() .* (R .- L)

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0 .= 1*x_1     # update value of all dimensions
                    break
                end
                # Value was not within slice, shrink the interval:
                #elseif all(x_1 .> x_0)
                if sign.(R .- x_0) == sign.(x_1 .- x_0)
                    R .= 1*x_1
                else
                    L .= 1*x_1
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        # Just overwrite the burnIn by using modular arithmetic:
        i = 1 + (ii-1)%N    # gives 1,2,3,...N,1,2,3,...,(N+burnIn)
        x_s[i,:] = x_1    # sample point, after one round of permutations
        pdflog_x_s[i] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations
    end
    return x_s, pdflog_x_s
end

# Diagonal Slice sampler using burn-ins, adapted from Jesper's code:
@everywhere function diagonalSliceSampling(pdf_log, x_0, w, N = 10_000, burnIn = 1_000; m = 100, verbose = true)
    # pre-burn-in run for getting a reasonable start point for the burn-in:
    x_s, l_pdf = mySliceSampler(pdf_log, x_0, w, burnIn; burnIn=0, m=m, verbose=verbose, msg="pre-burn-in")
    x_0 = x_s[argmax(l_pdf),:]   # Mode of pre-burn-in is new start point
    w = std(x_s, dims=1)  # a new guess for the window, my samples are column vectors

    println()
    # burn-in run for getting an idea about parameter covariation (bread-crums):
    x_s, l_pdf = mySliceSampler(pdf_log, x_0, w, burnIn; burnIn=0, m=m, verbose=verbose, msg="burn-in")

    println()
    # actual diagonal slice Sampling run, taking advantage of parameter Covariation
    C = cov(x_s)
    x_0 = x_s[argmax(l_pdf),:] # Mode of burn-in is new start point
    return mySliceSampler_Cov(pdf_log, x_0, C, N; burnIn=0, m=m, verbose=verbose, msg="MCMC run")
end

## Jespers HDI:
@everywhere function hdi(theta_samp,alpha=0.05)
    cred_mass = 1.0-alpha
    ci = zeros(2)
    if length(size(theta_samp))>1   # multidimensional theta
        K,N = size(theta_samp)
        cis = zeros(2,K)
        for k in 1:K    # for each dimension

            ts = theta_samp[k,:]
            sind = sortperm(ts)     # get series of indices that sorts theta values
            sts = ts[sind]  # sorted theta values

            # Shifting a 95% index block from far left to far right and use the
            # samllest interval that can be found during the shift
            N = length(sind)    # number of sample points
            length_ci = Inf     # start with inf HDI interval and shrink
            for i in 1:Int(floor(N*alpha))  # iterate over the 5% lowest values (low to high)
                i2 = Int(floor(N*cred_mass)+i)  # the 5% highest values (low to high), one 95% block away on the other side
                prop_ci = [sts[i],sts[i2]]      # proposed interval
                length_prop_ci = prop_ci[2]-prop_ci[1]
                if length_prop_ci < length_ci
                    ci = prop_ci
                    length_ci = ci[2]-ci[1]
                end
            end
            cis[:,k] = ci

        end
        return cis
    else
        N = length(theta_samp)

        ts = theta_samp
        sind = sortperm(ts)
        sts = ts[sind]

        N = length(sind)
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i],sts[i2]]
            length_prop_ci = prop_ci[2]-prop_ci[1]
            if length_prop_ci < length_ci
                ci = prop_ci
                length_ci = ci[2]-ci[1]
            end
        end
        return ci
    end
end

# A recursive application of the HDI algorithm to find the mode:
@everywhere function find_mode(sample_pts;runs=2)
    # re-apply the hdi algorithm, until there is only two points or so left
    # so the ci interval becomes smaller and smaller
    # in the end, return an average over the interval
    ci = zeros(2)
    i_left = 0
    i_right = 0

    # Make 50% credibility intervals
    alpha = 0.5;
    cred_mass = 1.0 - alpha;

    # sort the pts:
    sind = sortperm(sample_pts)
    sts = sample_pts[sind]
    N = length(sind)

    if N > runs
        # HDI algorithm, but return indices too:
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i], sts[i2]]
            length_prop_ci = prop_ci[2] - prop_ci[1]
            if length_prop_ci < length_ci
                i_left  = i
                i_right = i2
                ci = prop_ci
                length_ci = ci[2] - ci[1]
            end
        end

        # Recursion:
        ci = find_mode(sts[i_left:i_right],runs=runs)

    else
        ci = [sts[1],sts[end]]
    end

    return ci
end

# A function to make the plots with HDIs
function makeDistributionPlot(X, color="blue"; ann=true, offset=0.0, scale=1.0)
    #%% Mode, Mean, HDIs
    ω = mean(find_mode(X))
    μ_bar = mean(X);
    left,right = hdi(X);

    h = fit(Histogram, X, nbins=200)
    w = h.weights      # weights of each bar
    e = h.edges[1]     # edges of the bars (must be +1 more than bars)

    riemannSum = sum(w.*diff(e))
    w = w ./ riemannSum

    e2 = similar([e[:]; e[:]])  # doubling of the edge vector
    e2[2:2:end] = e             # even indices
    e2[1:2:end] = e             # odd indices
    e2 = e2[2:(end-1)]          # pop the ends

    w2 = similar([w[:]; w[:]])  # doubling of the weight vector
    w2[2:2:end] = w             # even indices
    w2[1:2:end] = w             # odd indices
    w2 = w2 .* scale .+ offset

    w_left = w[sum(e .< left)] * scale + offset  # comparison leaves at max left edge of bar that contains the limit
    w_right = w[sum(e .< right)] * scale + offset
    w_ω = w[sum(e .< ω)] * scale + offset
    w_μ = w[sum(e .< μ_bar)] * scale + offset

    #%% Plotting
    plot(e2, w2,
         linealpha=0.1, linecolor=color,
         fillrange=offset, fillalpha=0.4, fillcolor=color,
         legend=false)

    # or do it by histogram, then we also have vertical lines
    # specify statsplot?
    #histogram(X, bins=100, normalize=:pdf,
    #          color=color, alpha=0.4, linealpha=0.1, legend=false)   # Comes from StatsPlots now

    if ann
    plot!([(left,offset),(left,w_left)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (left, w_left, text("$(round(left,sigdigits=4))",color,rotation=90,:center,:left)))
    plot!([(right,offset),(right,w_right)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (right, w_right, text("$(round(right,sigdigits=4))",color,rotation=90,:center,:left)))
    plot!([(μ_bar,offset),(μ_bar,w_μ)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (μ_bar, 0.5*w_μ, text("mean $(round(μ_bar,sigdigits=4))",color,rotation=90,:bottom,:right)))
    plot!([(ω,offset),(ω,w_ω)],
        linewidth=1.5, color=color, linestyle=:dashdot,
        annotations = (ω, w_ω, text("mode $(round(ω,sigdigits=4))",color,rotation=90,:bottom,:right)))
    else
    plot!([(left,offset),(left,w_left)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(right,offset),(right,w_right)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(μ_bar,offset),(μ_bar,w_μ)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(ω,offset),(ω,w_ω)],
        linewidth=1.5, color=color, linestyle=:dashdot)
    end
end

function makeDistributionPlot!(X, color="blue"; ann=true, offset=0.0, scale=1.0)
    #%% Mode, Mean, HDIs
    ω = mean(find_mode(X))
    μ_bar = mean(X);
    left,right = hdi(X);

    h = fit(Histogram, X, nbins=200)
    w = h.weights      # weights of each bar
    e = h.edges[1]     # edges of the bars (must be +1 more than bars)

    riemannSum = sum(w.*diff(e))
    w = w ./ riemannSum

    e2 = similar([e[:]; e[:]])  # doubling of the edge vector
    e2[2:2:end] = e             # even indices
    e2[1:2:end] = e             # odd indices
    e2 = e2[2:(end-1)]          # pop the ends

    w2 = similar([w[:]; w[:]])  # doubling of the weight vector
    w2[2:2:end] = w             # even indices
    w2[1:2:end] = w             # odd indices
    w2 = w2 .* scale .+ offset

    w_left = w[sum(e .< left)] * scale + offset  # comparison leaves at max left edge of bar that contains the limit
    w_right = w[sum(e .< right)] * scale + offset
    w_ω = w[sum(e .< ω)] * scale + offset
    w_μ = w[sum(e .< μ_bar)] * scale + offset

    #%% Plotting
    plot!(e2, w2,
         linealpha=0.1, linecolor=color,
         fillrange=offset, fillalpha=0.4, fillcolor=color,
         legend=false)

    # or do it by histogram, then we also have vertical lines
    # specify statsplot?
    #histogram!(X, bins=100, normalize=:pdf,
    #          color=color, alpha=0.4, linealpha=0.1, legend=false)   # Comes from StatsPlots now

    if ann
    plot!([(left,offset),(left,w_left)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (left, w_left, text("$(round(left,sigdigits=4))",color,rotation=90,:center,:left)))
    plot!([(right,offset),(right,w_right)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (right, w_right, text("$(round(right,sigdigits=4))",color,rotation=90,:center,:left)))
    plot!([(μ_bar,offset),(μ_bar,w_μ)],
        linewidth=1.5, color=color, linestyle=:dash,
        annotations = (μ_bar, 0.5*w_μ, text("mean $(round(μ_bar,sigdigits=4))",color,rotation=90,:bottom,:right)))
    plot!([(ω,offset),(ω,w_ω)],
        linewidth=1.5, color=color, linestyle=:dashdot,
        annotations = (ω, w_ω, text("mode $(round(ω,sigdigits=4))",color,rotation=90,:bottom,:right)))
    else
    plot!([(left,offset),(left,w_left)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(right,offset),(right,w_right)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(μ_bar,offset),(μ_bar,w_μ)],
        linewidth=1.5, color=color, linestyle=:dash)
    plot!([(ω,offset),(ω,w_ω)],
        linewidth=1.5, color=color, linestyle=:dashdot)
    end
end

# Taken from Jesper:
# https://en.wikipedia.org/wiki/Correlation_coefficient
function ACF(x, k=1)
  x_meanshift = x .- mean(x)
  zxk = x_meanshift[k+1:end]
  zyk = x_meanshift[1:end-k]
  return sum(zxk.*zyk)/sqrt(sum(zxk.^2)*sum(zyk.^2))
end

# Taken from Jesper:
function acovlim(x;lim=0.05)
  k = 0
  rhos = []
  rho = 1
  while rho>lim
    rho = ACF(x,k)
    push!(rhos,rho)
    k += 1
  end
  return rhos
end

# Taken from Jesper:
# ess -- effective sample size (Kruschke 2014, page 184)
function ess(x)
    if typeof(x)==Vector{Float64}
        n = length(x)
        acf = acovlim(x)
        return n/(1+2*sum(acf[2:end]))
    else
        m,n = size(x)
        list = zeros(m)
        for i in 1:m
            acf = acovlim(x[i,:])
            list[i] = n/(1+2*sum(acf[2:end]))
        end
        return list
    end
end

################################################################################

#####
#%% Priors, Likelihoods, Posteriors ############################################
@everywhere function log_prior_flat(X::Float64)
    return 0.0  # or any other constant value, for the posterior it doesn't matter
end

@everywhere function log_prior_flat_positiveOnly(X::Float64)
    if X > 0
        return 0.0  # or any other constant value, for the posterior it doesn't matter
    else
        return -Inf
    end
    # Alternative of Jesper:
    #return log(tau > 0)
end

# Here we have to
# accomodate, that the means no longer are a constant, but an array, i.e.
# one mean for each theta
@everywhere function log_prior_theta(θ::Array{Float64},μ::Array{Float64},σ::Float64)
    # Log-normal distribution
    # Logarithmise the normal pdfs in series:
    # We can leave out the 1/sqrt(2π) since it does not change sampling from the posterior
    return sum(-log(σ) .- 0.5 * ((θ .- μ) ./ σ).^2)
end

@everywhere function log_lklhd_fct(J,jointParams::Array{Float64},zlogy::Array{Float64,1},zx::Array{Float64,1},ind::Array{Int64,1},K_j::Array{Int64,1})
    # return the joint lklhd of the given input data
    # Reparametrise, s.t. we sample eta ~ N(0,1) instead of theta. Theta is then reconstructed.
    η_0 = jointParams[1:J]
    η_1 = jointParams[(J+1):(2*J)]
    μ_0, ϕ_0, μ_1, ϕ_1, σ, τ_0, τ_1 = jointParams[(2*J+1):end]

    # Reconstruct thetas (span up by kids indicator K_j):
    θ_0 = (μ_0 .+ ϕ_0 * K_j) .+ τ_0 * η_0
    θ_1 = (μ_1 .+ ϕ_1 * K_j) .+ τ_1 * η_1

    jointLklhd = 0.0

    # the case of sigma <= 0 is already caught by the posterior main function

    # Logarithmise the normal pdfs in series:
    # summing up:
    # note, that we go through each observation, vectorising is hard:
    # zlogy ~ N(theta0 + theta1*zx, sigma)
    for i in 1:length(zlogy)
        j = ind[i]
        regMean = θ_0[j] + θ_1[j] * zx[i]   # regression mean
        jointLklhd += -log(σ) - 0.5 * ((zlogy[i] - regMean) / σ)^2
    end

    return jointLklhd
end

# Skipping the prior funciton speeds up the process.
# Still needs concretising by data
# Avoid using global vars, since slow
@everywhere function log_posterior_fct(J,jointParams::Array{Float64}, zlogy, zx, ind, child_j)
    # The joint log-pdf posterior of the model
    η_0 = jointParams[1:J]
    η_1 = jointParams[(J+1):(2*J)]
    μ_0, ϕ_0, μ_1, ϕ_1, σ, τ_0, τ_1 = jointParams[(2*J+1):end]

    #μ = μ_0 .+ ϕ .* K_j.*1.0    # current sampled μ for each individual -> this evaluation from A6 has now moved to the likelihood function due to reparametrisation!

    if (σ > 0) && (τ_0 > 0) && (τ_1 > 0)
        return log_lklhd_fct(J,jointParams::Array{Float64}, zlogy, zx, ind, child_j) +
               log_prior_theta(η_0,zeros(J),1.0) +  # eta0 ~ N(0,1)
               log_prior_theta(η_1,zeros(J),1.0) +  # eta1 ~ N(0,1)
               #log_prior_flat(μ_0) +   # is zero
               #log_prior_flat(μ_1) +   # is zero
               #log_prior_flat(ϕ_0) +   # is zero
               #log_prior_flat(ϕ_1) +   # is zero
               log_prior_flat_positiveOnly(σ) +
               log_prior_flat_positiveOnly(τ_0) +
               log_prior_flat_positiveOnly(τ_1)
    else
        return -Inf
    end
end

#%% Logarithmise and Transform Data ############################################
@everywhere J = maximum(ind)    # number of individuals
@everywhere I = length(y)       # number of observations

# Go to log-space:
@everywhere logy = log.(y)

@everywhere logMean = mean(logy)
@everywhere logStd = std(logy)

@everywhere trainMean = mean(x)
@everywhere trainStd = std(x)

# Mean-centre and scale each element:
@everywhere zlogy = (logy .- logMean) ./ logStd;

# Standardise training data:
@everywhere zx = (x .- trainMean) ./ trainStd;

#% Concretise the posterior by including the data #############################
@everywhere log_posterior(jointParams::Array{Float64}) =
            log_posterior_fct(J,jointParams::Array{Float64}, zlogy, zx, ind, child_j)

################################################################################
############################# MCMC #############################################
#%% Chain Setup
noParams = 2*J+7
noOfChains = 4
N = 4*10^6                # total number of samples
N_chain = convert(Int64, N / noOfChains)   # samples per chain
burnIn = 10^4   # the larger the burn-In the better the performance (better diagonal)
w = ones(Float64,noParams) * 0.1         # typical window size
m = 100                             # multiplier for maximum window size

#%% MCMC run, make several chains
chain = SharedArray{Float64}(N_chain,noParams,noOfChains);
chainPdf = SharedArray{Float64}(N_chain,noOfChains);

# Prerun to initialise the functions:
x_start = rand(Float64,noParams) * 0.5
@time trash1, trash2 = diagonalSliceSampling(log_posterior,x_start,w,100,1000,m=m,verbose=true);

# (using threads for parallel computing (other stuff didn't work)):
# Threads.@threads led to problems too
# now using @distributed, prefix with @sync to wait for completion
@sync @distributed for k in 1:noOfChains
    x_start = rand(Float64,noParams) * 0.5   # initial vector
    chain[:,:,k], chainPdf[:,k] = diagonalSliceSampling(log_posterior,x_start,w,N_chain,burnIn,m=m,verbose=true);
end

# My machine: 4 cores, i7-7500U 2,7 GHz x 2 each, on Linux
# N=10^5 takes around
# N=10^6 takes around 9-10 min
# N=4*10^6 takes around 40 min

################################################################################
############################# RESULTS ##########################################

##############
# Acess results and Transform back (undo mean-centering and scaling and go to non-log space)
# See derivation in PDF file instead

# Access the parameters and fuse the chains:
η_0 = Array{Float64,2}(undef,N,J)
η_1 = Array{Float64,2}(undef,N,J)
θ_0 = Array{Float64,2}(undef,N,J)
θ_1 = Array{Float64,2}(undef,N,J)
for j in 1:J
    η_0[:,j] = chain[:,j,:][:]
    η_1[:,j] = chain[:,J+j,:][:]
end

μ_0 = chain[:,2*J+1,:][:]
ϕ_0 = chain[:,2*J+2,:][:]
μ_1 = chain[:,2*J+3,:][:]
ϕ_1 = chain[:,2*J+4,:][:]
σ   = chain[:,2*J+5,:][:]
τ_0 = chain[:,2*J+6,:][:]
τ_1 = chain[:,2*J+7,:][:]

# This is slower than the for loop below:
# mu_intrcpt = repeat(μ_0,1,J) + repeat(ϕ_0,1,J) .* repeat(child_j,1,N)'
# mu_slp     = repeat(μ_1,1,J) + repeat(ϕ_1,1,J) .* repeat(child_j,1,N)'
# θ_0_t        = mu_intrcpt + η_0 .* repeat(τ_0,1,J)
# θ_1_t        = mu_slp     + η_1 .* repeat(τ_1,1,J)

for j = 1:J
    θ_0[:,j] = (μ_0 .+ ϕ_0 .* child_j[j]) .+ τ_0 .* η_0[:,j]
    θ_1[:,j] = (μ_1 .+ ϕ_1 .* child_j[j]) .+ τ_1 .* η_1[:,j]
end

### Exporting
# CSV.write(projDir*"/sliceSamples/export_samples_many.csv",
#          DataFrame([θ_0 θ_1 μ_0 ϕ_0 μ_1 ϕ_1 σ τ_0 τ_1]),
#          writeheader=true)

### Importing, if rerun
# read into a DataFrame:
# chain_csv = CSV.read(projDir*"/sliceSamples/export_samples_many.csv"; comment="#", normalizenames=true)
#
# θ_0 = Array{Float64,2}(undef,N,J)
# θ_1 = Array{Float64,2}(undef,N,J)
# for j in 1:J
#     θ_0[:,j] = chain_csv[:,j]
#     θ_1[:,j] = chain_csv[:,J+j]
# end
# μ_0 = chain_csv[:,2*J+1]
# ϕ_0 = chain_csv[:,2*J+2]
# μ_1 = chain_csv[:,2*J+3]
# ϕ_1 = chain_csv[:,2*J+4]
# σ   = chain_csv[:,2*J+5]
# τ_0 = chain_csv[:,2*J+6]
# τ_1 = chain_csv[:,2*J+7]
###

#####
# Autocorrelation and ESS
ess_θ_0 = Array{Float64,1}(undef,J)
ess_θ_1 = Array{Float64,1}(undef,J)
for j in 1:J
    ess_θ_0[j] = ess(θ_0[:,j])
    ess_θ_1[j] = ess(θ_1[:,j])
    println("ESS of θ_0_$j = ",ess_θ_0[j])
    println("ESS of θ_1_$j = ",ess_θ_1[j])
end
println("Minimal ESS of θ_0 is for individual $(argmin(ess_θ_0)) := ",ess_θ_0[argmin(ess_θ_0)])
println("Minimal ESS of θ_1 is for individual $(argmin(ess_θ_1)) := ",ess_θ_1[argmin(ess_θ_1)])
println("Maximal ESS of θ_0 is for individual $(argmax(ess_θ_0)) := ",ess_θ_0[argmax(ess_θ_0)])
println("Maximal ESS of θ_1 is for individual $(argmax(ess_θ_1)) := ",ess_θ_1[argmax(ess_θ_1)])

println("ESS of μ_0 = ",ess(μ_0))
println("ESS of ϕ_0 = ",ess(ϕ_0))
println("ESS of μ_1 = ",ess(μ_1))
println("ESS of ϕ_1 = ",ess(ϕ_1))
println("ESS of σ   = ",ess(σ))
println("ESS of τ_0 = ",ess(τ_0))
println("ESS of τ_1 = ",ess(τ_1))

#####
# Un-standardise:
# Individual level:
global θ_0_unscaled = logStd.*θ_0 .- logStd.*θ_1.*trainMean./trainStd .+ logMean
global θ_1_unscaled = logStd.*θ_1./trainStd

global σ_unscaled = σ .* logStd

# Group level (maybe not needed)
global μ_0_unscaled = logStd.*μ_0 .- logStd.*μ_1.*trainMean./trainStd .+ logMean
global ϕ_0_unscaled = logStd.*(ϕ_0 .- ϕ_1.*trainMean./trainStd)

global μ_1_unscaled = logStd.*μ_1./trainStd
global ϕ_1_unscaled = logStd.*ϕ_1./trainStd

global τ_0_unscaled = logStd.* sqrt.(τ_0.^2 .+ τ_1.^2 .* trainMean.^2 ./ trainStd.^2)
global τ_tilde_1_unscaled = sqrt.(2 .* trainMean) .* logStd .* τ_1 ./ trainStd
global τ_1_unscaled = logStd .* τ_1 ./ trainStd

#####
# Get into non-log space (generates function of input):
# Individual level:
function E_y_ind(x,j)
    # transform to zx:
    #zx = (x - trainMean) / trainStd
    unlogExpected = Array{Float64,1}(undef,N)
    unlogExpected = exp.(θ_0_unscaled[:,j] .+
                θ_1_unscaled[:,j] .* x .+
                0.5 .* σ_unscaled.^2);
    return unlogExpected
end

function E_y_group(x,isKid)
    # isKid is either 1 or 0
    # transform to zx:
    #zx = (x - trainMean) / trainStd
    unlogExpected = Array{Float64,1}(undef,N)
    unlogExpected = exp.(μ_0_unscaled + ϕ_0_unscaled .* isKid .+
                    (μ_1_unscaled .+ ϕ_1_unscaled .* isKid) .* x .+
                    0.5 .* (σ_unscaled.^2 .+ τ_0_unscaled.^2 .-
                            τ_tilde_1_unscaled.^2 .* x .+ τ_1_unscaled.^2 .* x^2))
    return unlogExpected
end

function curveSwarm(nrCurves,j,N,α,subpl=0)
    checkpoints = 100
    attempts = range(0,21,length=checkpoints)
    curves = Array{Float64,2}(undef,N,checkpoints)

    for (i, att) in enumerate(attempts)
        curves[:,i] = E_y_ind(att,j)
    end

    if subpl > 0
        plot!(legend=false,ylims=[100,800],subplot=subpl)
    else
        plot(legend=false,ylims=[100,800])
    end
    randIdces = Int.(ceil.(rand(nrCurves).*N))

    for i = 1:nrCurves
        if subpl > 0
            plot!(attempts,curves[randIdces[i],:],
                  color="blue",linealpha=α,subplot=subpl,
                  legend=false,ylims=[100,800])
        else
            plot!(attempts,curves[randIdces[i],:],
                  color="blue",linealpha=α,
                  legend=false,ylims=[100,800])
        end
    end
end

function curveSwarmGroup(nrCurves,isKid,N,α,subpl=0)
    checkpoints = 100
    attempts = range(0,21,length=checkpoints)
    curves = Array{Float64,2}(undef,N,checkpoints)

    for (i, att) in enumerate(attempts)
        curves[:,i] = E_y_group(att,isKid)
    end

    #if subpl > 0
    #    plot(legend=false,ylims=[100,800],subplot=subpl)
    #else
    #    plot(legend=false,ylims=[100,800])
    #end
    randIdces = Int.(ceil.(rand(nrCurves).*N))

    for i = 1:nrCurves
        if subpl > 0
            plot!(attempts,curves[randIdces[i],:],
                  color=RGB(isKid,0,0),linealpha=α,subplot=subpl,
                  legend=false,ylims=[100,800])
        else
            plot!(attempts,curves[randIdces[i],:],
                  color=RGB(isKid,0,0),linealpha=α,
                  legend=false,ylims=[100,800])
        end
    end
end

################################################################################
############################# TASKS ############################################

##### Task 3
for j in [1,3,4]
    ### Diagnostics to compare to Jesper:
    makeDistributionPlot(exp.(θ_1_unscaled[:,j]),"blue")
    plot!(xlabel="exp(theta_1[$j])")
    Plots.savefig(projDir*"/figs/th1$j-Slice.pdf")

    ### Expected reaction times:
    makeDistributionPlot(E_y_ind(1,j),"blue")
    makeDistributionPlot!(E_y_ind(5,j),"red")
    plot!(xlims=[200,800],grid=false,xlabel="E[y]-ind$j")
    Plots.savefig(projDir*"/figs/times-$j-Slice.pdf")

    ### Credible regression lines:
    curveSwarm(200,j,N,.15)
    col = RGB(child_j[j],0,0)
    scatter!(1:1:sum(ind.==j),y[ind.==j],
        markersize=10,markercolor=col,markerstrokecolor=col)
    plot!(grid=false,xlabel="attempt nr",ylabel="reaction time")
    Plots.savefig(projDir*"/figs/swarm-$j-Slice.pdf")
end

### All Regression lines:
plot(layout=(4, 9),size=(1500, 1000),
     legend=false,grid=false,xaxis=false,yaxis=false)
for j in 1:J
    curveSwarm(200,j,N,.15,j)
    if (j-1)%9==0
        myax=true
    else
        myax=false
    end
    col = RGB(child_j[j],0,0)
    scatter!(1:1:sum(ind.==j),y[ind.==j],
        markersize=4,markercolor=col,markerstrokecolor=col,
        subplot=j, ann=(12,700,"j = $j"),xaxis=true,yaxis=myax,
        framestyle=:axes)
end
Plots.savefig(projDir*"/figs/allcurves-Slice.pdf")

##### Diagnostics to compare to Jesper:
makeDistributionPlot(τ_0,"black")
Plots.savefig(projDir*"/figs/tau0-Slice.pdf")
makeDistributionPlot(τ_0_unscaled,"black")
Plots.savefig(projDir*"/figs/tau0-unsc-Slice.pdf")

makeDistributionPlot(σ,"black")
Plots.savefig(projDir*"/figs/sigma-Slice.pdf")
makeDistributionPlot(σ_unscaled,"black")
Plots.savefig(projDir*"/figs/sigma-unsc-Slice.pdf")

##### Swarm of the groups
attempts = range(0,22,length=200)
curve_mean_adult = mean.(E_y_group.(attempts,0))
curve_mean_kid = mean.(E_y_group.(attempts,1))
col = RGB.(child_i,0,0);

plot(attempts,curve_mean_adult,
    color=:black,ylims=[100,800],linewidth=8,linealpha=0.6)
plot!(attempts,curve_mean_kid,
      color=:red,ylims=[100,800],linewidth=8,linealpha=0.6)
curveSwarmGroup(200,0,N,.15)
curveSwarmGroup(200,1,N,.15)
plot!(grid=false,xlabel="attempt nr",ylabel="reaction time")
scatter!(x,y,markersize=2,markercolor=col,markerstrokecolor=col,alpha=0.2)
Plots.savefig(projDir*"/figs/swarm-groups-Slice.pdf")

### Chain diagnostics
# chstart = 650000
# chend = 655000
# plot(chstart:chend,θ_0[chstart:chend,:],legend=false,xlabel="sample nr",ylabel=L"\theta_{0_j}")
# Plots.savefig(projDir*"/figs/unstuckChain-theta0-Slice.pdf")
#
# plot(chstart:chend,θ_1[chstart:chend,:],legend=false,xlabel="sample nr",ylabel=L"\theta_{1_j}")
# Plots.savefig(projDir*"/figs/unstuckChain-theta1-Slice.pdf")
#
# plot(chstart:chend,σ[chstart:chend],legend=:topleft,xlabel="sample nr",label=L"\sigma")
# plot!(chstart:chend,ϕ_0[chstart:chend],xlabel="sample nr",label=L"\varphi_0")
# plot!(chstart:chend,ϕ_1[chstart:chend],xlabel="sample nr",label=L"\varphi_1")
# plot!(chstart:chend,μ_0[chstart:chend],xlabel="sample nr",label=L"\mu_0")
# plot!(chstart:chend,μ_1[chstart:chend],xlabel="sample nr",label=L"\mu_1")
# Plots.savefig(projDir*"/figs/unstuckChain-restPars-Slice.pdf")
# plot(chstart:chend,τ_0[chstart:chend],xlabel="sample nr",label=L"\tau_0")
# plot!(chstart:chend,τ_1[chstart:chend],xlabel="sample nr",label=L"\tau_1")
# Plots.savefig(projDir*"/figs/unstuckChain-tau-Slice.pdf")

makeDistributionPlot(τ_1,"red")
plot!(grid=false,xlabel=L"\tau_1")
Plots.savefig(projDir*"/figs/tau1-Slice.pdf")

# Convergence
chstart = 150000
chend = 160000
plot(grid=false)
for c = 0:noOfChains-1
    plot!(chstart:chend,τ_1[chstart+c*N_chain:chend+c*N_chain],
          xlabel="chain sample nr",
          label=latexstring("\\tau_1, chain_$(c+1)"),linealpha=0.4)
end
Plots.savefig(projDir*"/figs/convergence-Slice.pdf")

##### Useful posteriors
plot(grid=false,xlabel=L"\mu_1")
makeDistributionPlot!(μ_1,"black")
Plots.savefig(projDir*"/figs/mu1-Slice.pdf")
mean(exp.(μ_1_unscaled))

plot(grid=false,xlabel=L"\varphi_1")
makeDistributionPlot!(ϕ_1,"blue")
Plots.savefig(projDir*"/figs/phi1-Slice.pdf")
mean(exp.(ϕ_1_unscaled))

plot(grid=false,xlabel="Reaction time",xlims=(200, 600))
makeDistributionPlot!(E_y_group(1,0),"black",ann=false)
makeDistributionPlot!(E_y_group(1,1),"red",ann=false,scale=1,offset=0)
plot!(ann=(230,.005,"x = 1"))
makeDistributionPlot!(E_y_group(5,0),"black",ann=false,scale=1,offset=0.015)
makeDistributionPlot!(E_y_group(5,1),"red",ann=false,scale=1,offset=0.015)
plot!(ann=(230,.020,"x = 5"))
makeDistributionPlot!(E_y_group(10,0),"black",ann=false,scale=1,offset=0.03)
makeDistributionPlot!(E_y_group(10,1),"red",ann=false,scale=1,offset=0.03)
plot!(ann=(230,.035,"x = 10"))
Plots.savefig(projDir*"/figs/groupTimes-Slice.pdf")

##### Just for own reference and playing around with posterior predictions:
# zlogy_pred_kid(xin,Ξ) = mean(μ_0) .+ mean(ϕ_1) .+ (mean(μ_1) + mean(ϕ_1))*((xin - trainMean)/trainStd) .+
#                     sqrt(mean(σ)^2 + mean(τ_0)^2 + mean(τ_1)^2 * ((xin - trainMean)/trainStd)^2).*Ξ
#
# zlogy_pred_kid2(xin,Ξ,ix) = μ_0[ix] + ϕ_1[ix] + (μ_1[ix] + ϕ_1[ix])*((xin - trainMean)/trainStd) +
#                     sqrt(σ[ix]^2 + τ_0[ix]^2 + τ_1[ix]^2 * ((xin - trainMean)/trainStd)^2) * Ξ
#
# zlogy_pred_adult2(xin,Ξ,ix) = μ_0[ix] + μ_1[ix]*((xin - trainMean)/trainStd) +
#                     sqrt(σ[ix]^2 + τ_0[ix]^2 + τ_1[ix]^2 * ((xin - trainMean)/trainStd)^2) * Ξ
#
#
# atts = collect(1:21)
# plot()
# for i = 1:500
#     Ξ = randn()
#     randInd = max(Int(round(rand()*N)),1)
#     pKid = exp.((zlogy_pred_kid2.(atts,Ξ,randInd)) .* logStd .+ logMean)
#     pAdu = exp.((zlogy_pred_adult2.(atts,Ξ,randInd)) .* logStd .+ logMean)
#     plot!(atts,pKid,
#           color="red",linealpha=0.2,
#           legend=false)
#     plot!(atts,pAdu,
#           color="black",linealpha=0.2,
#           legend=false)
# end
#
# plot!()


text:
# This a dummy Julia module for demonstration purpose
module nimjlExample
  function customFunction(a::NamedTuple)
    println("From Julia -> args = ", a)
    # :nimTupKey is a trick to get the symbol name. It is equivalent to Symbol("nimTupKey1")
    nimTupKey1 = :nimTupKey1
    nimTupKey2 = :nimTupKey2

    # NamedTuple in Julia are immutable so it is necessary to create a new one
    result = (nimTupKey1 = a.nimTupKey1+1, nimTupKey2 = a.nimTupKey2+1)
    return result
  end

  mutable struct Foo
    x::Int
    y::Float64
    z::String
    # Nim initialize the Julia variable with empty constructor by default
    Foo() = new()
    Foo(x, y, z) = new(x, y, z)
  end

  function makeFoo()
    return Foo(1, -1.0, "Hello there")
  end

  function applyToFoo(foo::Foo)
    foo.x += 1
    foo.y *= 2/3
    foo.z *= " -> AppendMeBaby"
  end

  function returnDoubleValue()
    return (myint=36, mystr="Hello, world")
  end

  export returnDoubleValue
  export customFunction
  export Foo
  export makeFoo
  export applyToFoo

end

ext: jl

max_stars_repo_name: Clonkk/nimjl

max_issues_repo_name: Clonkk/nimjl

max_forks_repo_name: Clonkk/nimjl


text:
using Test

@testset "Dirichlet0 Boundary Conditions" begin
    N = 100
    h_inv = 1/(N-1)
    d_order = 2
    approx_order = 2
    x = collect(1:h_inv:N).^2
    A = DerivativeOperator{Float64}(2,2,1/99,length(x),:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    boundary_points = A.boundary_point_count
    res = A*x
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 2.0*ones(N - sum(boundary_points));

    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N).^4 .- 2*collect(1:1.0:N).^3 .+ collect(1:1.0:N).^2;
    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    boundary_points = A.boundary_point_count
    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-1; # Float64 is less stable

    A = DerivativeOperator{BigFloat}(d_order,approx_order,one(BigFloat),N,:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    y = convert(Array{BigFloat, 1}, y)
    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^(-approx_order);
end

@testset "Dirichlet0 Boundary Conditions irregular grid" begin
    # tests as for DerivativeOperator, but with arbitrary grid (increasing step size):
    N = 1000
    d_order = 2
    approx_order = 2
    x = (1.:1.:N)*2.0 #erf.(linspace(-2,2,N))*N/2;
    dx = diff(x)
    y = x.^2

    A = FiniteDifference{Float64}(2,2,dx,length(x),:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    boundary_points = A.boundary_point_count
    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 2.0*ones(N - sum(boundary_points));

    # N = 1000
    d_order = 4
    approx_order = 10
    # dx = rand(N)
    # x = [0.0; cumsum(dx)]
    y = @. x.^4 - 2*x.^3 + x.^2
    A = FiniteDifference{Float64}(d_order,approx_order,dx,N,:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    boundary_points = A.boundary_point_count
    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^(-1); # Float64 is less stable

    A = FiniteDifference{BigFloat}(d_order,approx_order,BigFloat.(dx),N,:Dirichlet0,:Dirichlet0,BC=(0.0,0.0))
    y = convert(Array{BigFloat, 1}, y)
    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^(-approx_order);

end


text:
facts("Exchange responses for backtesting") do
  context("Software-implemented limit order") do
    plim = 150.0
    buyord = TradingLogic.Order("buysome", 10, plim,
                                :buy, :limit, :pending, "aux")
    sellord = TradingLogic.Order("sellsome", 10, plim,
                                 :sell, :limit, :pending, "aux")
    @fact TradingLogic.plimitcheck(buyord, 175.0) --> false
    @fact TradingLogic.query_orderstatus(buyord, 175.0) --> :pending
    @fact TradingLogic.plimitcheck(buyord, 125.0) --> true
    @fact TradingLogic.query_orderstatus(buyord, 125.0) --> :complete
    @fact TradingLogic.plimitcheck(sellord, 175.0) --> true
    @fact TradingLogic.query_orderstatus(sellord, 175.0) --> :complete
    @fact TradingLogic.plimitcheck(sellord, 125.0) --> false
    @fact TradingLogic.query_orderstatus(sellord, 125.0) --> :pending
  end
  context("Market order and order submission") do
    ordid = "markord"
    mord = TradingLogic.Order(ordid, 10, NaN,
                              :buy, :market, :pending, "aux")
    # buy no matter how high
    @fact TradingLogic.query_orderstatus(mord, 1e7) --> :complete
    # simulated order submission
    @fact TradingLogic.submit_ordernew(mord, false) --> "FAIL"
    @fact TradingLogic.submit_ordernew(mord, true) --> ordid
  end
  context("Cancel pending limit-order") do
    limord = TradingLogic.Order("buysome", 10, 100.0,
                                :buy, :limit, :pending, "aux")
    @fact TradingLogic.submit_ordercancel(limord) --> true
    # can not cancel non-pending order
    @fact TradingLogic.submit_ordercancel(
      TradingLogic.emptyorder()) --> false
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Pyramid(p1, p2, p3, p4, p5)

A pyramid with points `p1`, `p2`, `p3`, `p4`, `p5`.
"""
struct Pyramid{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polyhedron{Dim,T}
  vertices::V
end


text:
function mass_matrix_filter_entries!(
    mass_matrix,
    filt::F,
    global_index::Base.ImmutableDict{Symbol, Int64},
) where {F <: PSY.Filter}
    @debug "Using default mass matrix entries $F"
    return
end

function mdl_filter_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    dynamic_device::DynamicWrapper{PSY.DynamicInverter{C, O, IC, DC, P, PSY.LCLFilter}},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}

    #external_ix = get_input_port_ix(dynamic_device, PSY.LCLFilter)
    #θ_oc = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[Vr_inv_var]
    V_tI = inner_vars[Vi_inv_var]
    Vr_cnv = inner_vars[Vr_cnv_var]
    Vi_cnv = inner_vars[Vi_cnv_var]

    #Get parameters
    filter = PSY.get_filter(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)
    ωb = 2 * pi * f0
    lf = PSY.get_lf(filter)
    rf = PSY.get_rf(filter)
    cf = PSY.get_cf(filter)
    lg = PSY.get_lg(filter)
    rg = PSY.get_rg(filter)
    basepower = PSY.get_base_power(dynamic_device)
    sys_Sbase = get_system_base_power(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.LCLFilter)

    #Define internal states for filter
    internal_states = @view device_states[local_ix]
    Ir_cnv = internal_states[1]
    Ii_cnv = internal_states[2]
    Vr_filter = internal_states[3]
    Vi_filter = internal_states[4]
    Ir_filter = internal_states[5]
    Ii_filter = internal_states[6]

    #Inputs (control signals) - N/A

    #Compute 6 states ODEs (D'Arco EPSR122 Model)
    #Inverter Output Inductor (internal state)
    #𝜕id_c/𝜕t
    output_ode[local_ix[1]] = (
        ωb / lf * Vr_cnv - ωb / lf * Vr_filter - ωb * rf / lf * Ir_cnv +
        ωb * ω_sys * Ii_cnv
    )
    #𝜕iq_c/𝜕t
    output_ode[local_ix[2]] = (
        ωb / lf * Vi_cnv - ωb / lf * Vi_filter - ωb * rf / lf * Ii_cnv -
        ωb * ω_sys * Ir_cnv
    )
    #LCL Capacitor (internal state)
    #𝜕vd_o/𝜕t
    output_ode[local_ix[3]] =
        (ωb / cf * Ir_cnv - ωb / cf * Ir_filter + ωb * ω_sys * Vi_filter)
    #𝜕vq_o/𝜕t
    output_ode[local_ix[4]] =
        (ωb / cf * Ii_cnv - ωb / cf * Ii_filter - ωb * ω_sys * Vr_filter)
    #Grid Inductance (internal state)
    #𝜕id_o/𝜕t
    output_ode[local_ix[5]] = (
        ωb / lg * Vr_filter - ωb / lg * V_tR - ωb * rg / lg * Ir_filter +
        ωb * ω_sys * Ii_filter
    )
    #𝜕iq_o/𝜕t
    output_ode[local_ix[6]] = (
        ωb / lg * Vi_filter - ωb / lg * V_tI - ωb * rg / lg * Ii_filter -
        ωb * ω_sys * Ir_filter
    )

    #Update inner_vars
    inner_vars[Vr_filter_var] = Vr_filter
    inner_vars[Vi_filter_var] = Vi_filter

    #Compute current from the inverter to the grid
    I_RI = (basepower / sys_Sbase) * [Ir_filter; Ii_filter]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]
    return
end

function mdl_filter_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ω_sys::ACCEPTED_REAL_TYPES,
    dynamic_device::DynamicWrapper{PSY.DynamicInverter{C, O, IC, DC, P, PSY.RLFilter}},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}
    #Obtain inner variables for component
    basepower = PSY.get_base_power(dynamic_device)
    sys_Sbase = get_system_base_power(dynamic_device)
    ratio_power = basepower / sys_Sbase

    #Obtain parameters
    filt = PSY.get_filter(dynamic_device)
    rf = PSY.get_rf(filt)
    lf = PSY.get_lf(filt)

    Vr_cnv = inner_vars[Vr_cnv_var]
    Vi_cnv = inner_vars[Vi_cnv_var]
    Vr_inv = inner_vars[Vr_inv_var]
    Vi_inv = inner_vars[Vi_inv_var]

    #Compute output currents
    if lf != 0.0 || rf != 0.0
        Zmag_squared = rf^2 + lf^2
        Ir_filt = (1.0 / Zmag_squared) * ((Vr_cnv - Vr_inv) * rf + (Vi_cnv - Vi_inv) * lf)
        Ii_filt = (1.0 / Zmag_squared) * ((Vi_cnv - Vi_inv) * rf - (Vr_cnv - Vr_inv) * lf)
    else
        #Obtain converter
        IS.@assert_op Vr_cnv == Vr_inv
        IS.@assert_op Vi_cnv == Vi_inv
        converter = PSY.get_converter(dynamic_device)
        R_source = PSY.get_R_source(converter)
        X_source = PSY.get_X_source(converter)
        Z_source_sq = R_source^2 + X_source^2
        I_aux_r = (R_source * Vr_cnv + X_source * Vi_cnv) / Z_source_sq
        I_aux_i = (R_source * Vi_cnv - X_source * Vr_cnv) / Z_source_sq
        Ir_filt = inner_vars[Ir_cnv_var] - I_aux_r
        Ii_filt = inner_vars[Ii_cnv_var] - I_aux_i
    end
    #Update Inner Vars
    inner_vars[Ir_inv_var] = Ir_filt
    inner_vars[Ii_inv_var] = Ii_filt

    #Update current
    current_r[1] += ratio_power * Ir_filt
    current_i[1] += ratio_power * Ii_filt
    return
end


text:
# TODO: add comments (including assumptions, and quirks of subtype)

struct SimplifiedMethod
	name::String
	sig::DataType
end

datatype_name(dt::DataType) = dt.name

primary_datatype(dt::DataType) = datatype_name(dt).primary

is_abstract(dt::DataType) = dt.abstract

function isa_tuple_adjust(x::Any, t::Type)
	if !isa(t, DataType) || (t.name != Tuple.name) || (t === Tuple)
		isa(x, t)
	else
		if !isa(x, Tuple)
			return false
		end
		xidx = 1
		pidx = 1
		ps = collect(t.parameters)
		while xidx <= length(x)
			if pidx > length(ps)
				return false
			end
			if !isa(x[xidx], isa(ps[pidx], TypeVar) ? ps[pidx].ub : (ps[pidx].name == Vararg.name) ? ps[pidx].parameters[1] : ps[pidx])
				return false
			end
			xidx += 1
			if ps[pidx].name != Vararg.name
				pidx += 1
			end
		end
		(pidx == length(ps)+1) || ((pidx <= length(ps)) && (ps[pidx].name == Vararg.name))
	end
end

function datatype_tree_ascend(dt::DataType, dttree::Dict{DataType, Vector{DataType}} = Dict{DataType, Vector{DataType}}())
	if dt != Any
		superdt = supertype(dt)
		v = get!(dttree, superdt, Vector{DataType}[])
		if !(dt in v)
			push!(v, dt)
		end
		datatype_tree_ascend(superdt::DataType, dttree)
	end
	dttree
end

function datatype_tree_descend(dt::DataType, dttree::Dict{DataType, Vector{DataType}} = Dict{DataType, Vector{DataType}}(); subtypefn=subtypes)
	v = get!(dttree, dt, Vector{DataType}[])
	for subdt in subtypefn(dt)
		if subdt != Any
			if !(subdt in v)
				push!(v, subdt)
			end
			datatype_tree_descend(subdt::DataType, dttree; subtypefn=subtypefn)
		end
	end
	dttree
end

function datatype_tree(dt::DataType, dttree::Dict{DataType, Vector{DataType}} = Dict{DataType, Vector{DataType}}(); subtypefn=subtypes)
	datatype_tree_ascend(dt, dttree)
	datatype_tree_descend(dt, dttree; subtypefn=subtypefn)
	dttree
end

function datatype_tree(dts::Vector{DataType}, dttree::Dict{DataType, Vector{DataType}} = Dict{DataType, Vector{DataType}}(); subtypefn=subtypes)
	for dt in dts
		datatype_tree(dt::DataType, dttree; subtypefn=subtypefn)
	end
	dttree
end

function nonabstract_descendents(dt::DataType, descendents::Vector{DataType} = Vector{DataType}(); subtypefn=subtypes)
	if !is_abstract(dt)
		if !(dt in descendents)
			push!(descendents, dt)
		end
	else
		for subdt in subtypefn(dt)
			if subdt != Any
				nonabstract_descendents(subdt::DataType, descendents; subtypefn=subtypefn)
			end
		end
	end
	descendents
end

function nonabstract_descendents(dts::Vector{DataType}, descendents::Vector{DataType} = Vector{DataType}(); subtypefn=subtypes)
	for dt in dts
		nonabstract_descendents(dt::DataType, descendents; subtypefn=subtypefn)
	end
	descendents
end

function extract_primary_datatypes(dt::DataType, pdts::Vector{DataType} = Vector{DataType}())
	pdt = primary_datatype(dt)
	if !(pdt in pdts)
		push!(pdts, pdt)
	end
	for p in dt.parameters
		extract_primary_datatypes(p, pdts)
	end
	pdts
end

function extract_primary_datatypes(ts::Vector{Type}, pdts::Vector{DataType} = Vector{DataType}())
	for t in ts
		extract_primary_datatypes(t, pdts)
	end
	pdts
end

function extract_primary_datatypes(u::Union, pdts::Vector{DataType} = Vector{DataType}())
	for t in u.types
		extract_primary_datatypes(t, pdts)
	end
	pdts
end

function extract_primary_datatypes(tv::TypeVar, pdts::Vector{DataType} = Vector{DataType}())
	extract_primary_datatypes(tv.ub, pdts)
	extract_primary_datatypes(tv.lb, pdts) # since Union{} is a Union rather than a DataType, and it contains no types, it will contribute nothing
	pdts
end

function extract_primary_datatypes(t::Any, pdts::Vector{DataType} = Vector{DataType}())
	# do nothing (catch all needed for e.g. array dimensions which are numeric values)
	pdts
end

function merge_datatypes_up(dts::Vector{DataType})
	if isempty(dts)
		return dts
	else
		if length(dts) > 1
			for i in 2:length(dts)
				if dts[1] <: dts[i]
					return merge_datatypes_up(dts[2:end])
				end
				if dts[i] <: dts[1]
					dts[i] = dts[1]
					return merge_datatypes_up(dts[2:end])
				end
			end
		end
		[dts[1]; merge_datatypes_up(dts[2:end])]
	end
end

function merge_datatypes_down(dts::Vector{DataType})
	if isempty(dts)
		return dts
	else
		if length(dts) > 1
			for i in 2:length(dts)
				if dts[1] <: dts[i]
					dts[i] = dts[1]
					return merge_datatypes_down(dts[2:end])
				end
				if dts[i] <: dts[1]
					return merge_datatypes_down(dts[2:end])
				end
			end
		end
		[dts[1]; merge_datatypes_down(dts[2:end])]
	end
end

is_partially_supported(dt::DataType, supporteddts::Vector{DataType}) = begin
	any(supporteddt->typeintersect(supporteddt, dt) !== Union{}, supporteddts) &&
	all(p -> is_partially_supported(p, supporteddts), dt.parameters)
end

is_partially_supported(u::Union, supporteddts::Vector{DataType}) = all(p -> is_partially_supported(p, supporteddts), u.types)
	# TODO really could do any, but for the purposes of the DataGenerator it will pick a Union at random, so really need all to be supported

is_partially_supported(tv::TypeVar, supporteddts::Vector{DataType}) = is_partially_supported(tv.ub, supporteddts)
	# TODO lb?

# Switchoff becouse TypeConstructor not recomended in julia after 1.0
#is_partially_supported(tc::TypeConstructor, supporteddts::Vector{DataType}) = is_partially_supported(tc.body, supporteddts)

is_partially_supported(x::Any, supporteddts::Vector{DataType}) = true # for type parameters such as integers etc.

# is_partially_supported(ts::Vector{Type}, supporteddts::Vector{DataType}) = all(t -> is_partially_supported(t, supporteddts), ts)

function is_partially_supported(m::SimplifiedMethod, supporteddts::Vector{DataType})
	@assert m.sig <: Tuple
	all(p->is_partially_supported(p, [supporteddts; Vararg]), m.sig.parameters)
end

partially_supported_constructor_methods(dt::DataType, supporteddts::Vector{DataType}) =
	filter(m->is_partially_supported(m, supporteddts), simplified_methods(dt)) # TODO restrict by module?  # TODO other filters (e.g. deprecated)

is_fully_supported(dt::DataType, supporteddts::Vector{DataType}) = begin
	(((dt == Any) && !isempty(supporteddts)) # Any is implicitly supported if there is at least datatype
		|| any(supporteddt->typeintersect(supporteddt, dt) === dt, supporteddts)) && # NB may be bug with ==: Union{} == dt for all dt
	all(p -> is_fully_supported(p, supporteddts), dt.parameters)
end

is_fully_supported(u::Union, supporteddts::Vector{DataType}) = all(p -> is_fully_supported(p, supporteddts), u.types)
	# TODO really could do any, but for the purposes of the DataGenerator it will pick a Union at random, so really need all to be supported

is_fully_supported(tv::TypeVar, supporteddts::Vector{DataType}) = is_fully_supported(tv.ub, supporteddts)
	# TODO lb?

# Switchoff becouse TypeConstructor not recomended in julia after 1.0
#is_fully_supported(tc::TypeConstructor, supporteddts::Vector{DataType}) = is_fully_supported(tc.body, supporteddts)

is_fully_supported(x::Any, supporteddts::Vector{DataType}) = true # for type parameters such as integers etc.

# is_fully_supported(ts::Vector{Type}, supporteddts::Vector{DataType}) = all(t -> is_fully_supported(t, supporteddts), ts)

function is_fully_supported(m::SimplifiedMethod, supporteddts::Vector{DataType})
	@assert m.sig <: Tuple
	all(p->is_fully_supported(p, [supporteddts; Vararg]), m.sig.parameters)
end

fully_supported_constructor_methods(dt::DataType, supporteddts::Vector{DataType}) =
	filter(m->is_fully_supported(m, supporteddts), simplified_methods(dt)) # TODO restrict by module?  # TODO other filters (e.g. deprecated)

type_as_parseable_string(tv::TypeVar) = "TypeVar(Symbol(\"" * string(tv.name) * "\")," * type_as_parseable_string(tv.lb) * "," * type_as_parseable_string(tv.ub) * "," * (tv.bound ? "true" : "false") * ")"

type_as_parseable_string(dt::DataType) = string(dt.name) * (dt === primary_datatype(dt) ? "" : "{" * join(map(p -> type_as_parseable_string(p), dt.parameters), ",") * "}") # note: === rather than == on testing primary datatypes is necessary to ensure parameters match (e.g. bound/unbound TypeVars)

# Switchoff becouse TypeConstructor not recomended in julia after 1.0
#type_as_parseable_string(tc::TypeConstructor) = type_as_parseable_string(tc.body)

type_as_parseable_string(u::Union) = "Union" * "{" * join(map(p -> type_as_parseable_string(p), u.types), ",") * "}" # note: Union{} is distinct from Union, so always output the curly braces

type_as_parseable_string(x::Any) = string(x)

replace_typevars_with_ub(tv::TypeVar) = tv.ub

replace_typevars_with_ub(dt::DataType) = primary_datatype(dt){map(p -> replace_typevars_with_ub(p), dt.parameters)...}

replace_typevars_with_ub(u::Union) = Union{map(p -> replace_typevars_with_ub(p), u.types)...}

replace_typevars_with_ub(x::Any) = x

simplified_methods(dt::DataType) = map(m->simplify_method(m, dt), methods(dt))

function simplify_method(m::Method, dt::DataType)
	if m.name == :Type
		@assert m.sig <: Tuple
		methodname = type_as_parseable_string(dt)
		# if isa(m.tvars, TypeVar)
		# 	if m.tvars.ub != Any
		# 		methodname *= "{" * type_as_parseable_string(m.tvars.ub) * "}"
		# 	end
		# else
		# 	for tvar in m.tvars
		# 		if tvar.ub != Any
		# 			methodname *= "{" * type_as_parseable_string(tvar.ub) * "}"
		# 		end
		# 	end
		# end
		SimplifiedMethod(methodname, Tuple{m.sig.parameters[2:end]...})
		# println("For method $(m): converting to $(sm.name) with sig $(sm.sig)")
	else
		SimplifiedMethod(string(m.name), m.sig)
	end
end


text:
# ---
# title: 258. Add Digits
# id: problem258
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/add-digits/description/>
# hidden: true
# ---
#
# Given a non-negative integer `num`, repeatedly add all its digits until the
# result has only one digit.
#
# **Example:**
#
#
#
#     Input: 38
#     Output: 2
#     Explanation: The process is like: 3 + 8 = 11, 1 + 1 = 2.
#                  Since 2 has only one digit, return it.
#
#
# **Follow up:**
# Could you do it without any loop/recursion in O(1) runtime?
#
#
## @lc code=start
using LeetCode

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


text:
import LazySets.use_precise_ρ

"""
    DiscretePost

Abstract supertype of all discrete post operators.

### Notes

All discrete post operators should provide the following method, in addition
to those provided for general post operators:
```julia
tube⋂inv(𝒫::DiscretePost, reach_tube::Vector{<:AbstractReachSet{<:LazySet, N}},
         invariant, start_interval) where {N}
```
"""
abstract type DiscretePost <: PostOperator end

function postprocess(𝒫,
                     HS,
                     post_jump,
                     options,
                     waiting_list,
                     passed_list,
                     target_loc_id,
                     jumps)
    fixpoint_strategy = options[:fixpoint_check]

    if fixpoint_strategy == :eager
        # eager fixpoint checking
        post_jump_filtered =
            filter(x -> !isfixpoint(𝒫, x, passed_list, target_loc_id),
                   post_jump)
    else
        post_jump_filtered = post_jump
    end

    if (isempty(post_jump_filtered))
        # fixpoint found or transition can never be taken
        return
    end

    # apply clustering
    clustered = cluster(𝒫, post_jump_filtered, options)

    # push new sets after jump (unless a fixpoint is detected)
    for reach_set in clustered
        if fixpoint_strategy != :none
            if fixpoint_strategy == :lazy &&
                    isfixpoint(𝒫, reach_set, passed_list, target_loc_id)
                continue
            end
            push!(passed_list[target_loc_id], reach_set)
        end
        push!(waiting_list, (target_loc_id, reach_set, jumps))
    end
end

function cluster(𝒫::DiscretePost,
                 reach_sets::Vector{RSN},
                 options::Options) where {SN, RSN<:AbstractReachSet{SN}}
    clustering_strategy = options[:clustering]
    oa = 𝒫.options[:overapproximation]
    if clustering_strategy == :none
        # no clustering, keeping original set
        return reach_sets
    elseif clustering_strategy == :none_oa
        # no clustering but overapproximation
        return [RSN(overapproximate(set(reach_set), oa),
        time_start(reach_set), time_end(reach_set)) for reach_set in reach_sets]
    elseif clustering_strategy == :chull
        # cluster all sets in a convex hull and overapproximate that set
        chull = ConvexHullArray([set(reach_set) for reach_set in reach_sets])
        chull_oa = overapproximate(chull, oa)
        return [RSN(chull_oa, time_start(reach_sets[1]),
                    time_end(reach_sets[end]))]
    end
end

function isfixpoint(𝒫::DiscretePost,
                    reach_set::RSN,
                    passed_list,
                    loc_id
                   ) where {SN, RSN<:AbstractReachSet{SN}}
    @assert passed_list != nothing
    if isassigned(passed_list, loc_id)
        for other_reach_set in passed_list[loc_id]
            if set(reach_set) ⊆ set(other_reach_set)
                info("found a fixpoint in some reach tube")
                return true
            end
        end
        return false
    else
        passed_list[loc_id] = Vector{RSN}()
        return false
    end
end

# default: always apply line search
function use_precise_ρ(𝒫::DiscretePost,
                             cap::Intersection{N})::Bool where N<:Real
    return true
end

# --- default methods for handling assignments ---

# default implementation: use 'apply' from MathematicalSystems
function apply_assignment(𝒫::DiscretePost,
                          constrained_map::AbstractMap,
                          R⋂G::LazySet;
                          kwargs...)
    return apply(constrained_map, R⋂G)
end

# for reset maps: return a lazy ResetMap from LazySets
function apply_assignment(𝒫::DiscretePost,
                          constrained_map::ConstrainedResetMap,
                          R⋂G::LazySet;
                          kwargs...)
    return LazySets.ResetMap(R⋂G, constrained_map.dict)
end


text:
using Printf

configuration = [[ 30; 110.0; 9700.0; 14.0;11.0; -9.0], # 1
                 [ 40;  50.0; 2600.0;  5.0; 1.0; -7.0], # 2
                 [ 50;  80.0; 4900.0; 20.0; 2.0; -5.5], # 3
                 [ 60;  30.0; 8600.0; 30.0; 2.0; -5.0], # 4
                 [ 70; 150.0; 5000.0;  7.0; 2.0; -9.0], # 5
                 [ 80;  90.0; 3600.0;  4.0;20.0; -7.5], # 6
                 [160;  90.0; 3600.0;  4.0;20.0; -7.5], # 7
                 [120;  90.0; 3600.0;  4.0;20.0; -7.5], # 8
                 [180;  70.0; 9000.0; 10.0;10.0; -5.5], # 9
                 [100;  50.0; 7000.0; 12.0;12.0; -7.0], # 10
                 [200;  50.0; 7000.0; 12.0;12.0; -6.0], # 11
                 [150;  56.0; 2200.0; 12.0;12.0; -5.0], # 12
                 [100;  50.0; 4000.0;  3.0; 3.0; -7.0], # 13
                 [100;  50.0; 5000.0;  3.0; 3.0; -6.0], # 14
                 [100;  50.0; 3000.0;  2.0; 4.0; -6.0], # 15
                 [100;  50.0; 5500.0;  2.0; 5.0; -5.0], # 16
                 [300;  30.0; 6600.0; 30.0; 2.0; -7.0], # 17
                 [100;  30.0; 6700.0; 15.0; 2.0; -6.0], # 18
                 [100;  70.0; 6700.0; 11.0; 2.0; -8.0], # 19
                 [100;  70.0; 6700.0; 11.0; 2.0;-10.0], # 20
                 [100;  70.0; 6700.0; 16.0; 2.0;-11.0], # 21
                 [110; 110.0; 6000.0; 19.0; 2.0; -9.0], # 22
                 [100;  80.0; 7100.0;  5.0; 2.0; -8.0]] # 23

for i in 1:length(configuration)

    N  = round(Int, configuration[i][1])
    SN = configuration[i][2]
    T  = configuration[i][3]
    F  = configuration[i][4]
    S  = configuration[i][5]
    LS = configuration[i][6]

    @show configuration[i]

    t, flux, σ, idx, regflux, regσ, line, _ = simulatedata(N=N, SN=SN, Temperature=T, offset=F, scale=S, logs=LS, seed=i);

    figure(3)
    savefig(@sprintf("GP_total_fluxes_%d.png", i))
    figure(4)
    savefig(@sprintf("GP_flux_plot_%d.png", i))

    #posterior = bpca(regflux', (regσ)',maxiter=600)
    posterior = gpmodelVI(t + 0*randn(size(t)), flux, σ.^2, idx; logs=LS, maxiter=2_000)

    g = unitgalaxyvector(T)*S

    randx = noisyintersectionvi(posterior=posterior, g=g)

    xavg = mean([randx() for i in 1:10_000])

    figure(20)
    cla()

    subplot(611)
    cla()
    plt.hist([randx()[1] for i=1:100_000], 30)
    plot(g[1], 0, "ro", label="true coordinate 1")
    plot(xavg[1], 0, "ko", label="xavg 1")
    legend()

    subplot(612)
    cla()
    plt.hist([randx()[2] for i=1:100_000], 30)
    plot(g[2], 0, "ro", label="true coordinate 2")
    plot(xavg[2], 0, "ko", label="xavg 2")
    legend()

    subplot(613)
    cla()
    plt.hist([randx()[3] for i=1:100_000], 30)
    plot(g[3], 0, "ro", label="true coordinate 3")
    plot(xavg[3], 0, "ko", label="xavg 3")
    legend()

    subplot(614)
    cla()
    plt.hist([randx()[4] for i=1:100_000], 30)
    plot(g[4], 0, "ro", label="true coordinate 4")
    plot(xavg[4], 0, "ko", label="xavg 3")
    legend()

    subplot(615)
    cla()
    plt.hist([randx()[5] for i=1:100_000], 30)
    plot(g[5], 0, "ro", label="true coordinate 5")
    plot(xavg[5], 0, "ko", label="xavg 3")
    legend()

    subplot(616)
    cla()
    plt.hist([randx()[6] for i=1:100_000], 30)
    plot(g[6], 0, "ro", label="true coordinate 6")
    plot(xavg[6], 0, "ko", label="xavg 3")
    legend()

    savefig(@sprintf("GP_hist_%d.png", i))

end


text:
using GPUArrays
using Base.Test, GPUArrays.TestSuite

# It's kind of annoying to make FillArrays only a test dependency
# so for texting the conversion to GPUArrays of shaped iterators,
# I just copied the core types from FillArrays:s

abstract type AbstractFill{T, N} <: AbstractArray{T, N} end
@inline function Base.getindex(F::AbstractFill, k::Integer)
    @boundscheck checkbounds(F, k)
    getindex_value(F)
end
@inline function Base.getindex(F::AbstractFill{T, N}, kj::Vararg{<:Integer, N}) where {T, N}
    @boundscheck checkbounds(F, kj...)
    getindex_value(F)
end
Base.IndexStyle(F::AbstractFill) = IndexLinear()
struct Fill{T, N} <: AbstractFill{T, N}
    value::T
    size::NTuple{N, Int}
end
getindex_value(x::Fill) = x.value
@inline Base.size(F::Fill) = F.size
struct Eye{T} <: AbstractMatrix{T}
    size::NTuple{2, Int}
end
Base.size(E::Eye) = E.size
@inline function Base.getindex(E::Eye{T}, k::Integer, j::Integer) where T
    @boundscheck checkbounds(E, k, j)
    ifelse(k == j, one(T), zero(T))
end

function run_construction(Typ)
    @testset "Construction" begin
        constructors(Typ)
        conversion(Typ)
        value_constructor(Typ)
        iterator_constructors(Typ)
    end
end

function constructors(Typ)
    @testset "similar + constructor" begin
        for T in supported_eltypes()
            B = Typ{T}(10)
            @test size(B) == (10,)
            @test eltype(B) == T

            B = Typ{T}(10, 10)
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = Typ{T}((10, 10))
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = similar(B, Int32, (11, 15))
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(B, T)
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(B, (5,))
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(B, 7)
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(Typ{Int32}, (11, 15))
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(Typ{Int32, 2}, T, (11, 15))
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(Typ{T}, (5,))
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(Typ{T}, 7)
            @test size(B) == (7,)
            @test eltype(B) == T
        end
    end
end

function conversion(Typ)
    @testset "conversion" begin
        for T in supported_eltypes()
            Bc = round.(rand(10, 10) .* 10.0)
            B = Typ{T}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10)
            B = Typ(Bc)
            @test size(B) == (10,)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10, 10)
            B = Typ{T, 2}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(Int32, 3, 3, 3)
            B = convert(Typ{T, 3}, Bc)
            @test size(B) == (3, 3, 3)
            @test eltype(B) == T
            @test Array(B) ≈ Bc
        end
    end
end

function value_constructor(Typ)
    @testset "value constructor" begin
        for T in supported_eltypes()
            x = zeros(T, 2, 2)
            x1 = zeros(Typ{T}, 2, 2)
            x2 = zeros(Typ{T}, (2, 2))
            x3 = zeros(Typ{T, 2}, (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            x = ones(T, 2, 2)
            x1 = ones(Typ{T}, 2, 2)
            x2 = ones(Typ{T}, (2, 2))
            x3 = ones(Typ{T, 2}, (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            fill!(x, 0)

            x1 = fill(Typ{T}, 0f0, 2, 2)
            x2 = fill(Typ{T}, T(0), (2, 2))
            x3 = fill(Typ{T, 2}, T(0), (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            fill!(x1, 2f0)
            x2 = fill!(Typ{Int32}((4, 4, 4)), 77f0)
            @test all(x-> x == 2f0, Array(x1))
            @test all(x-> x == Int32(77), Array(x2))

            x = eye(T, 2, 2)

            x1 = eye(Typ{T, 2}, 2, 2)
            x2 = eye(Typ{T}, (2, 2))
            x3 = eye(Typ{T, 2}, (2, 2))

            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x
        end
    end
end
function iterator_constructors(Typ)
    @testset "iterator constructors" begin
        for T in supported_eltypes()
            @test Typ(Fill(T(0), (10,))) == zeros(Typ{T}, 10)
            @test Typ(Fill(T(0), (10, 10))) == zeros(Typ{T}, 10, 10)
            if T <: Real
                x = Typ{Float32}(Fill(T(0), (10, 10)))
                @test eltype(x) == Float32
                @test Typ(Eye{T}((10, 10))) == eye(Typ{T}, 10, 10)
                x = Typ{Float32}(Eye{T}((10, 10)))
                @test eltype(x) == Float32
            end
        end
    end
end


text:
## dataframes.jl DataFrames support for ROC
## (c) 2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

using DataFrames: AbstractDataFrame, DataFrame

## We also implement most TNT arguments here.

## Accept a DataFrame with :target and :score columns to most functions
"""
TNT(::DataFrame; optional-args) extracts target and non-target scores from a dataframe.  The optional arhuments encode which colums are used for what purpose. `score` (default `:score`) is the name of the column containing the classifier's scores.  `target` (default `:target`), of type `Bool`, determines whether this is a target score or not.
"""
function TNT(x::AbstractDataFrame; score=:score, target=:target)
    ok = .!(ismissing.(x[!, target]) .| ismissing.(x[!, score]))
    t = BitArray(x[!, target][ok])
    s = collect(skipmissing(x[!, score][ok]))
    TNT(s[t], s[.!t])
end

for f in (:eer, :eerch, :auc, :cllr, :mincllr)
    @eval ($f)(tnt::TNT) = ($f)(tnt.tar, tnt.non)
    @eval ($f)(x::AbstractDataFrame; score=:score, target=:target) = ($f)(TNT(x, score=score, target=target))
end

roc(x::AbstractDataFrame; score=:score, target=:target, laplace=false, collapse=true) = roc(TNT(x, score=score, target=target), laplace=laplace, collapse=collapse)
roc(tnt::TNT; laplace=false, collapse=true) = roc(tnt.tar, tnt.non, laplace=laplace, collapse=collapse)

for f in (:ber, :minber)
    @eval ($f)(tnt::TNT, plo) = ($f)(tnt.tar, tnt.non, plo)
    @eval ($f)(x::AbstractDataFrame, plo; score=:score, target=:target) = ($f)(TNT(x, score=score, target=target), plo)
end

dcf(x::AbstractDataFrame; score=:score, target=:target, d::DCF=getdcf(), thres=-plo(d), norm=false) = dcf(TNT(x, score=score, target=target), d=d, thres=thres, norm=norm)

mindcf(x::AbstractDataFrame; score=:score, target=:target, d::DCF=getdcf(), norm=false) = mindcf(TNT(x, score=score, target=target), d=d, norm=norm)

import DataFrames.DataFrame
"""
`Dataframe(::Roc)` converts a `Roc` object in a datframe with columns:

 - `pfa` the false alarm rate

 - `pmiss` the miss rate

 - `thres` the thereshold, separating this line's pfa and pmiss from the next

 - `chull` indicating if this point is on the convex hull of the ROC curve

 - `llr` the optimal log-likelihood-ratio score for all data points contributing to the ROC line segment from this line to the next
"""
DataFrame(r::Roc) = DataFrame(pfa=r.pfa, pmiss=r.pmiss, thres=[Array(r.θ); missing], chull=r.ch, llr=[Array(r.llr); missing])


text:
function calculate_pi(NBIN::Int= 1000000000)::Float64
    step = 1/NBIN
    sum = 0
    for i=1:NBIN
        x = (i+0.5)*step;
        sum += 4 / (1 + x^2)
    end
    pi = sum*step
    return pi
end

if length(ARGS)>0
    @time @show calculate_pi(parse(Int, ARGS[1]))
else
    @time @show calculate_pi()
end

max_stars_repo_name: naezzell/accelqat

max_issues_repo_name: naezzell/accelqat

max_forks_repo_name: naezzell/accelqat


text:
# pandigital number
p = [1,2,3,4,5,6,7,8,9]

function next_pandigital(p::Array{Int,1})
    N = 9
    i = N-1
    while i > 0
        # find trigger p[end-i]
        trigger = p[i]
        if trigger < p[i+1]
            s = p[i:end]    # subarray of numbers that have to be swapped
            a = minimum(s[s .> trigger])   # smallest number that is larger than p[end-i]
            sort!(filter!(x->x≠a,s))

            # do the swapping
            p[i] = a
            p[i+1:end] = s
            break
        else
            i -= 1
        end

        if i == 0   # pandigital number is 987654321
            p = p[end:-1:1]
        end
    end
    return p
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
# ENV["JULIA_DEBUG"] = Main
struct MyToken
	token::String
	tipo::Symbol
end

struct MyExpr
	call::Symbol
	args::Vector{Any}
end

abstract type MyExprAST end
struct NumberExprAST <: MyExprAST
	val::Float64
end

struct HeadFunExprAST <: MyExprAST
	name::Symbol
	args::Vector{Symbol}
end

struct IfExprAST <: MyExprAST
    cond::MyExprAST
    then::MyExprAST
    elsee::MyExprAST
end

struct BodyFunExprAST <: MyExprAST
	head::HeadFunExprAST
	body::MyExprAST
end

struct BodyTailExprAST <: MyExprAST
	head::HeadFunExprAST
	body::IfExprAST
end
BodyTailExprAST(expr::BodyFunExprAST) = BodyTailExprAST(expr.head, expr.body)

struct DefVarExprAST <: MyExprAST
	varname::Symbol
	init::MyExprAST
end

struct VarExprAST <: MyExprAST
	name::Symbol
end

struct ReduceExprAST <: MyExprAST
	fun::Symbol
    args::Vector{MyExprAST}
end

struct CallExprAST <: MyExprAST
    call::Symbol
    args::Vector{MyExprAST}
end

struct BoolOpExprAST <: MyExprAST
	op::Symbol
	rhs::MyExprAST
	lhs::MyExprAST
end

const cursor = Ref{Int}()

const specials = ['!', '@', '#', '$', '^', '&', '*', '-', '_', '+', '=', '<', '>', '/', '?', ';']

myisletter(char) = char in 'a':'z' || char in specials
myisnumber(char) = char in '0':'9'
myisacii(char) = myisletter(char)||myisnumber(char)
myisatom(token) = token.tipo == :name || token.tipo == :number || token.tipo == :string
myisopen(token) = token.tipo == :open
myisclosed(token) = token.tipo == :closed
myisname(token) = token.tipo == :name
myistring(token) = token.tipo == :string
myisnumber(token::MyToken) = token.tipo == :number
function lexer(str)
	if str[1] == '"'
		ind = findfirst('"',str[2:end])
		if ind === nothing
			throw(ErrorException("String não foi fechada"))
		end
		return MyToken(str[1:ind], :string)
	end
	if str[1] == '('
		return MyToken("(", :open)
	end
	if str[1] == ')'
		return MyToken(")", :closed)
	end
	if str[1] == '-' || str[1] == '+'
		if str[2] == ' '
			return MyToken(string(str[1]), :name)
		end
		ind = findfirst(!myisnumber,str[2:end]) + 1
		if ind === nothing
			error("+ and - can't be used in the start of identifiers")
		end
		if str[ind] == '.'
			ind2 = findfirst(!myisnumber,str[ind+1:end])
			if ind === nothing
				return MyToken(str, :number)
			end
			return MyToken(str[1:ind+ind2-1], :number)
		end
		return MyToken(str[1:ind-1], :number)
		ind
		if all(myisnumber.(str[2:end]))
			return MyToken(str, :number)
		else
			error("+ and - can't be used in the start of identifiers")
		end
	end
	if myisnumber(str[1])
		ind = findfirst(!myisnumber,str)
		if ind === nothing
			return MyToken(str, :number)
		end
		if str[ind] == '.'
			ind2 = findfirst(!myisnumber,str[ind+1:end])
			if ind === nothing
				return MyToken(str, :number)
			end
			return MyToken(str[1:ind+ind2-1], :number)
		end
		return MyToken(str[1:ind-1], :number)
	end
	if myisletter(str[1])
		ind = findfirst(!myisacii,str)
		if ind === nothing
			return MyToken(str, :name)
		end
		return MyToken(str[1:ind-1], :name)
	end
	throw(ErrorException("Invalid char $(str[1])"))
end

function increment(token)
	cursor[] += length(token.token)
end

function gettoken(str)
	token = lexer(str[cursor[]:end])
	increment(token)
	return token
end

function walkspace(str)
	ind = findfirst(x-> x != ' '&& x != '\n',str[cursor[]:end])
	if ind === nothing
		return nothing
	end
	cursor[] += ind -1
end

function pusharg!(expr::MyExpr, arg)
	push!(expr.args, arg)
	nothing
end

function parse_token(token::MyToken)
	if myisname(token)
		return Symbol(token.token)
	elseif myisnumber(token)
		return parse(Float64,token.token)
	elseif myistring(token)
		# return token.token
		throw(ArgumentError("Strings não implementadas"))
	else
		throw(ArgumentError("Erro ao processar token"))
	end
end

function expr_machine(str)
	walkspace(str)
	token = gettoken(str)
	@debug "Token gerado" token
	if !myisname(token)
		throw(ErrorException("Syntax error $(token.token)"))
	end
	expr = MyExpr(parse_token(token), [])
	@debug "expr: Estado 1 para 2"
	walkspace(str)
	token2 = gettoken(str)
	@debug "Token gerado" token2
	if myisopen(token2)
		@debug "expr: Estado 2 para 3 e chamou submáquina de expressões"
		pusharg!(expr, expr_machine(str))
		@debug "expr: A submáquina de expressões terminou e foi para o estado 4"
	elseif myisatom(token2)
		@debug "expr: Estado 2 para 4"
		pusharg!(expr, parse_token(token2))
	else
		throw(ErrorException("Syntax error $(token2.token)"))
	end
	walkspace(str)
	token3 = gettoken(str)
	@debug "Token gerado" token3
	while !myisclosed(token3)
		if myisopen(token3)
			@debug "expr: Estado 4 para 3 e chamou submáquina de expressões"
			pusharg!(expr, expr_machine(str))
			@debug "expr: A submáquina de expressões terminou e foi para o estado 4"
		elseif myisatom(token3)
			@debug "expr: Estado 4 para 4"
			pusharg!(expr, parse_token(token3))
		else
			throw(ErrorException("Syntax error $(token3.token)"))
		end
		walkspace(str)
		token3 = gettoken(str)
		@debug "Token gerado" token3
	end
	return expr
end

function root_machine(str)
	walkspace(str)
	token = gettoken(str)
	@debug "Token gerado" token
	if myisopen(token)
		@debug "root: Chamou a submáquina de expressões"
		expr = expr_machine(str)
		@debug "root: A submáquina de expressões terminou"
	return expr
	end
	if myisatom(token)
		@debug "root: Estado 1 pra 2 com atómo"
		return parse_token(token)
	end

	throw(ErrorException("Syntax error $(token.token)"))
end

function global_machine(str)
	@debug "Inicio do processo"
	@debug "Código a ser processado" str
	cursor[] = 1
	if length(str) == 1
		return [root_machine(str)]
	end

	exprs = Any[]
	while cursor[] < length(str)
		@debug "Chamou a máquina da raiz da expressão"
		push!(exprs, root_machine(str))
		@debug "Retornou da máquina da raiz da expressão"
	end
	return exprs
end


text:
# ------------ Call sites ------------ #

@inline function (ctx::ScoreContext)(c::typeof(markov),
                                     addr::Address,
                                     call::Function,
                                     len::Int,
                                     args...)
    visit!(ctx, addr => 1)
    ss = get_subselection(ctx, addr)
    ret, w = score(get_sub(ss, 1), call, args...)
    v_ret = Vector{typeof(ret)}(undef, len)
    v_ret[1] = ret
    increment!(ctx, w)
    for i in 2:len
        visit!(ctx, addr => i)
        ret, w = score(get_sub(ss, i), call, v_ret[i-1]...)
        v_ret[i] = ret
        increment!(ctx, w)
    end
    return v_ret
end


text:
using Oceananigans.Grids: topology, min_Δx, min_Δy, min_Δz

"Returns the time-scale for advection on a regular grid across a single grid cell."
function cell_advection_timescale(u, v, w, grid)

    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)

    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)

    return min(Δx/umax, Δy/vmax, Δz/wmax)
end

cell_advection_timescale(model) =
    cell_advection_timescale(model.velocities.u.data.parent,
                             model.velocities.v.data.parent,
                             model.velocities.w.data.parent,
                             model.grid)


text:
include("pick_inN.jl")
include("pick_inS.jl")
include("driversBracket.jl")

# Version inspired by JCG algo 6.6 (Fletcher-Lemaréchal) which avoids computing
# derivatives when not needed, e.g. if pick_in is the basic mid point.

export bracket

"""  `bracket(ϕ, a, b, stp, α, β, pick_in, best)`
   - ϕ is a linemodel
   - minimization in the interval [a,b], default [0,∞].
      - [a,b] is an initial interval quaranteed to contain a local minimum; if not satisfied, the function exits. In particular, ϕ'(a) must be negative.
   - stp : stopping object; atol and rtol are bypassed and the conditions below apply
   - α<0<β is the (possibly asymmetric) stopping tolerance, default β=-α=1e-8
   - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<h(a)
   - pick_in : function to pick a point in an interval;
      -  several variants are proposed, default enhanced secant
   - best option to bracket the best candidate or any minimizer such that h(t*)<h(a), default = true

`bracket` is a composite implementation of
   - a bracketing scheme
   - an interval reduction scheme using a pick_in function to select a point in a given interval
   - enhanced pick_in variants involving a polynomial interpolation (Newton like)
"""
function bracket(ϕ       :: OneDModel;
                 a       :: T = 0.0 ,   b  :: T = Inf,
                 #stp     :: AbstractStopping = NLPStopping(ϕ, OneDAtX(a)),
                 stp     :: AbstractStopping =
                 NLPStopping(ϕ, OneDAtX(a, zeros(2)),
                             tol_check=(a,b,c) -> ones(2)),
                 α       :: T = -1e-8,  β  :: T = 1e-8,
                 pick_in :: Function = pick_ins2,  # defaults to secant's enhancement
                 best    :: Bool = true,           # aim for the lowest minimizer
                 #kwargs...
                 ) where T <: AbstractFloat

    ϕa, dϕa = NaN,  NaN
    if isa(ϕ, LSModel) # then ϕ.f₀ and ϕ.g₀ are initialized
        @assert a==0
        ϕa, dϕa = T(0), (1-ϕ.τ₀)*ϕ.g₀
    else  # require a function and derivative computation
        ϕa, dϕa = obj(ϕ, a), grad(ϕ, a)
    end

    tref = a
    ϕref = ϕa
    dϕref = dϕa

    ϕ₀ = ϕa
    ####################
    # à revoir, mais ça semble OK, on tente de trouver un meilleur minimiseur dans l'algo,
    # mais le critère d'arrêt ne teste que la descente par rapport à t=a.
    # on stoppe donc lorsque ϕ(t)<ϕ(a)   et   α ≤ dϕ(t) ≤  β
    #
    # Ce serait plus propre que ϕa=0 toujours...autrement, avec les erreurs d'arrondi,
    # le truc de tricher la valeur de ϕ courante ne marche possiblement pas... voir(*) plus bas
    ####################
    stp.meta.optimality_check = (p::OneDModel, s::AbstractState) -> [s.fx; s.gx] #
    stp.meta.atol = T(0)
    stp.meta.rtol = T(0)
    stp.meta.tol_check_neg = (a,b,c) -> [-Inf; α]
    stp.meta.tol_check     = (a,b,c) -> [  ϕ₀; β]

    ϕb, dϕb = NaN, NaN;    if isfinite(b)  ϕb, dϕb = obj(ϕ, b), grad(ϕ, b)  end

    tₘ, ϕₘ, dϕₘ = a, ϕa, dϕa # best objective so far

    OK, status = Check_coherent(a, b, ϕa, dϕa, ϕb, dϕb, α, β)
    if !OK
        @warn string("in one D minimization: ",status)
        @error "in one D minimization: non coherent"
        stp.meta.suboptimal = true
        return stp
    end

    OK = update_and_start!(stp, x = a, fx=ϕa, gx = dϕa)

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"h(t)", :dual=>"h'(t)", :step=>"a", :slope=>"b"))

    @debug  "α = ", α, "   β = ", β
    dϕt = Inf
    forward = true
    initInf!()

    while !OK
        t = pick_in(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        # avoid computing derivatives as long as descent is not obtained
        ϕt = obj(ϕ, t); dϕt = Inf # grad(ϕ, t)

        # for as cheap a linesearch as possible, avoid to compare with ϕₘ, keep on
        #                                                 comparing to 0 (Armijo condition)
        if best          # rather compare with ϕₘ to aim the lowest minimizer possible.
            tref = tₘ
            ϕref = ϕₘ
            dϕref = dϕₘ
        end

        descent = ϕt < ϕref
        if ~descent # double check for eventual roundoff error dominating the descent computation
            if (abs(ϕref - ϕt) < 1e-7) # heuristic, apply at roughly half the machine precision
                dϕt = grad(ϕ, t)       # not sure that the grad is already available
                val_descent = (dϕref+dϕt)*(t - tref)/2
                    descent = val_descent < 0
            end
            # if descent, cheat and set the value of ϕt to its "descent" value       (*)
            if descent ϕt = ϕref + val_descent end
            @debug "best  ", ϕ₀, ϕₘ, ϕt, tₘ, t, dϕₘ, dϕt, descent
        end

        if descent
            # t and ϕt becomes the best guess so far
            tₘ, ϕₘ = t, ϕt
            dϕt = grad(ϕ, t)
            @debug " New grad dϕt = ", dϕt, "  t = ", t, "α = ", α, "   β = ", β
            @debug "           ϕt = ", ϕt, "  ϕ₀ = ", ϕ₀, " ϕₘ = ", ϕₘ
            dϕₘ = dϕt
            if dϕt > β
                b, ϕb, dϕb = t, ϕt, dϕt
                forward = false
            elseif dϕt < α
                a, ϕa, dϕa = t, ϕt, dϕt
                forward = true
            else  # should stop
                @debug "   stop  dϕt =  $dϕt "
            end
        elseif forward
            b, ϕb = t, ϕt
        else  #backward
            a, ϕa = t, ϕt
        end
        @debug " avant stop   ϕt = $ϕt    ϕₘ = $ϕₘ, dϕₘ = $dϕₘ , dϕt = $dϕt "
        OK = update_and_stop!(stp, x=tₘ, fx=ϕₘ, gx=dϕₘ)
        @info log_row(Any[stp.meta.nb_of_stop, ϕt, dϕt, a, b])

    end
    if ϕₘ > ϕ₀ @warn "Increase in line minimization" end
    return stp
end

# Tricky implementation to have increasing jumps when b is infinite
let f_mult_inc = 5.0, inc = 1.0
    global function setInf!(v)
        inc = v
    end

    global function initInf!(;val :: Float64 = 1.0)
        setInf!(val)
        resetN3!() # resets the history for quintic Newton interpolation
        resetSec!()# resets the history for secant variants
    end

    global function getInf()
        return f_mult_inc, inc
    end
end

"""
    pick_inInf is the common increase function when b=Inf
"""
function pick_inInf(a)
    f_mult_inc, inc = getInf()
    anew = a + inc
    inc *= f_mult_inc
    setInf!(inc)
    return anew
end

"""
pick_inB is the basic picking scheme which picks the middle point in [a, b].

"""
function pick_inB(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)

    if b < Inf
        return (a + b)/2.0
    else
        return pick_inInf(a)
    end
end

""" Check that the given interval [a,b] is guaranteed to contain a solution
    such that α ≤ ϕ' ≤ β"""
function Check_coherent(a, b, ϕa, dϕa, ϕb, dϕb, α, β)

    OK = true
    status = :OK
    #  Check that the given interval [a,b] is guaranteed to contain a solution
    if b <= a
        OK = false
        status = :BadInterval
    elseif (dϕa >= 0.0)
        OK = false
        status = :NonDecrease
        #elseif !implies(b<Inf, (dϕb > β) || (ϕb >= ϕa))
        #    OK = false
        #    status = :NoSolutionGaranty
    elseif b<Inf
        if ϕb<ϕa
            if dϕb < β
                if dϕb > α
                    OK = true
                    status = :solved  #  b satisfies   (ϕb<ϕa)  &&  (α < dϕb < β)
                elseif dϕb < 0
                    OK = false
                    status = :NoSolutionWaranty
                end
            end
        end
    end

    return OK, status
end


text:
using Test, Random, LinearAlgebra, SparseArrays
using Yao
using YaoExtensions
using QuAlgorithmZoo

@testset "PhaseEstimation" begin
    include("PhaseEstimation.jl")
end

@testset "hamiltonian solvers" begin
    include("hamiltonian_solvers.jl")
end

@testset "hadamard test" begin
    include("HadamardTest.jl")
end