text:
function problem_generator_mw(
                            nprob::Int64,
                            n::Int64,
                            m::Int64,
                            rsp::Bool;
                            unconstrained::Bool=false,
                            x::Vector{Float64}=( NaN * ones(Float64, n) )
                            )

    @assert 0 ≤ nprob ≤ 22 "The number of the problem must satisfy 0 ≤ nprob ≤ 22."
    @assert length(x) == n "The dimension of vetor x must be n."

    # Initializes some useful variables and constants
    c_lower = - 1.0e20
    c_upper = 1.0e20

    if isnan( x[1] )

        usual_init_point = true

    end

    if nprob == 1

        # Linear function - Full Rank.

        @assert m ≥ n "The dimensions of the problem must satisfy m ≥ n."

        if usual_init_point

            @. x = 1.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = 0.75 * ones(Float64, n)
            u = 15.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)
        c = 2.0 / m

        for i = 1:n

            fmin[i] = x -> ( x[i] - c * sum(x) - 1.0 ) ^ 2.0

        end

        for i = (n+1):m

            fmin[i] = x -> ( - c * sum(x) - 1.0 ) ^ 2.0

        end

    elseif nprob == 2

        # Linear function - Rank 1.

        @assert m ≥ n "The dimensions of the problem must satisfy m ≥ n."

        if usual_init_point

            @. x = 1.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - ones(Float64, n)
            u = 10.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        function f2(x)

            s = 0.0

            for i=1:length(x)

                s += i * x[i]

            end

            return s

        end

        for i=1:m

            fmin[i] = x -> ( i * f2(x) - 1.0 ) ^ 2.0

        end

    elseif nprob == 3

        # Linear function - Rank 1 with zero columns and rows.

        @assert m ≥ n "The dimensions of the problem must satisfy m ≥ n."

        if usual_init_point

            @. x = 1.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = 2.0 * ones(Float64, n)
            u = 7.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        function f3(x)

            s = 0.0

            for i=2:(length(x)-1)

                s += i * x[i]

            end

            return s

        end

        fmin[1] = x -> 1.0
        fmin[m] = x -> 1.0
        for i = 2:(m-1)

            fmin[i] = x -> ( ( i - 1.0 ) * f3(x) - 1.0 ) ^ 2.0

        end

    elseif nprob == 4

        # Rosenbrock function.

        @assert ( n == 2 ) && ( m == 2 ) "The dimensions of the problem must satisfy n = m = 2."

        if usual_init_point

            @. x = [-1.2, 1.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 5.0 * ones(Float64, n)
            u = zeros(Float64, n)

        end

        fmin = Array{Function}(undef, m)
        fmin[1] = x -> ( 10.0 * ( x[2] - x[1] ^ 2.0 ) ) ^ 2.0
        fmin[2] = x -> ( 1.0 - x[1] ) ^ 2.0

    elseif nprob == 5

        # Helical Valley function

        @assert ( n == 3 ) && ( m == 3 ) "The dimensions of the problem must satisfy n = m = 3."

        if usual_init_point

            @. x = [-1.0, 0.0, 0.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 10.0 * ones(Float64, n)
            u = - ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        function f5(x)

            if x[1] > 0.0

                return atan( x[2] / x[1] ) / ( 2.0 * pi )

            elseif x[1] < 0.0

                return atan( x[2] / x[1] ) / ( 2.0 * pi ) + 0.5

            else

                return 0.25

            end

        end

        fmin[1] = x -> ( 10.0 * ( x[3] - 10.0 * f5(x) ) ) ^ 2.0
        fmin[2] = x -> ( 10.0 * ( sqrt( x[1] ^ 2.0 + x[2] ^ 2.0 ) - 1.0 ) ) ^ 2.0
        fmin[3] = x -> ( x[3] ) ^ 2.0

    elseif nprob == 6

        # Powell Singular function

        @assert ( n == 4 ) && ( m == 4 ) "The dimensions of the problem must satisfy n = m = 4."

        if usual_init_point

            @. x = [3.0, -1.0, 0.0, 1.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - ones(Float64, n)
            u = 10.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)
        fmin[1] = x -> ( x[1] + 10.0 * x[2] ) ^ 2.0
        fmin[2] = x -> 5.0 * ( x[3] - x[4] ) ^ 2.0
        fmin[3] = x -> ( x[2] - 2.0 * x[3] ) ^ 4.0
        fmin[4] = x -> 10.0 * ( x[1] - x[4] ) ^ 4.0

    elseif nprob == 7

        # Freudenstein and Roth function

        @assert ( n == 2 ) && ( m == 2 ) "The dimensions of the problem must satisfy n = m = 2."

        if usual_init_point

            @. x = [0.5, -2.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 3.0 * ones(Float64, n)
            u = 3.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)
        fmin[1] = x -> ( - 13.0 + x[1] + ( ( 5.0 - x[2] ) * x[2] - 2.0 ) * x[2] ) ^ 2.0
        fmin[2] = x -> ( - 29.0 + x[1] + ( ( x[2] + 1.0 ) * x[2] - 14.0 ) * x[2] ) ^ 2.0

    elseif nprob == 8

        # Bard function

        @assert ( n == 3 ) && ( m == 15 ) "The dimensions of the problem must satisfy n = 3 e m = 15."

        if usual_init_point

            @. x = 1.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = 0.5 * ones(Float64, n)
            u = 10.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        y = [
                0.14, 0.18, 0.22, 0.25, 0.29,
                0.32, 0.35, 0.39, 0.37, 0.58,
                0.73, 0.96, 1.34, 2.10, 4.39
                ]

        for i=1:m

            v = 16.0 - i
            w = min( i, v )

            fmin[i] = x -> ( y[i] - ( x[1] + ( i ) / ( v * x[2] + w * x[3] ) ) ) ^ 2.0

        end

    elseif nprob == 9

        # Kowalik and Osborne function

        @assert ( n == 4 ) && ( m == 11 ) "The dimensions of the problem must satisfy n = 4 e m = 11."

        if usual_init_point

            @. x = [0.25, 0.39, 0.415, 0.39]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = [0.0, 0.0, 0.415, 0.0]
            u = [0.275, 0.4, 1.0, 0.4]

        end

        fmin = Array{Function}(undef, m)
        y1 = [ 0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246 ]
        y2 = [ 4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625 ]

        for i=1:m

            fmin[i] = x -> ( y1[i] - ( x[1] * ( y2[i] ^ 2.0 + y2[i] * x[2] ) ) / ( y2[i] ^ 2.0 + y2[i] * x[3] + x[4] ) ) ^ 2.0

        end

    elseif nprob == 10

        # Meyer function

        @assert ( n == 3 ) && ( m == 16 ) "The dimensions of the problem must satisfy n = 3 e m = 16."

        if usual_init_point

            @. x = [0.02, 4000.0, 250.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 100.0 * ones(Float64, n)
            u = [100.0, 4000.0, 1000.0]

        end

        fmin = Array{Function}(undef, m)
        y = [
                34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0,
                8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0 ]

        for i=1:m

            fmin[i] = x -> ( x[1] * exp( x[2] / ( 45.0 + 5.0 * i + x[3] ) ) - y[i] ) ^ 2.0

        end

    elseif nprob == 11

        # Watson function

        @assert ( 2 ≤ n ≤ 31 ) && ( m == 31 ) "The dimensions of the problem must satisfy 2 ≤ n ≤ 31 e m = 31."

        if usual_init_point

            @. x = 0.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = ones(Float64, n)
            u = 10.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        function f11a(x, ti)

            s = 0.0

            for j = 2:length(x)

                s += ( j - 1.0 ) * x[j] * ti ^ ( j - 2.0 )

            end

            return s

        end

        function f11b(x, ti)

            s = 0.0

            for j = 1:length(x)

                s += x[j] * ti ^ ( j - 1.0 )

            end

            return s

        end

        for i = 1:29

            ti = i / 29.0

            fmin[i] = x -> ( f11a(x, ti) - ( f11b(x, ti) ) ^ 2.0 - 1.0 ) ^ 2.0

        end

        fmin[30] = x -> ( x[1] ) ^ 2.0
        fmin[31] = x -> ( x[2] - x[1] ^ 2.0 - 1.0 ) ^ 2.0

    elseif nprob == 12

        # Box three-dimensional function

        @assert ( n == 3 ) && ( m ≥ n ) "The dimensions of the problem must satisfy n = 3 e m ≥ n."

        if usual_init_point

            @. x = [0.0, 10.0, 20.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = [-1.0, -50.0, 15.0]
            u = [50.0, 10.0, 50.0]

        end

        fmin = Array{Function}(undef, m)

        for i=1:m

            ti = 0.1 * i

            fmin[i] = x -> ( exp( - ti * x[1] ) - exp( - ti * x[2] ) - x[3] * ( exp( - ti ) - exp( - 10.0 * ti ) ) ) ^ 2.0

        end

    elseif nprob == 13

        # Jennrich and Sampson function

        @assert ( n == 2 ) && ( m ≥ n ) "The dimensions of the problem must satisfy n = 2 e m ≥ n."

        if usual_init_point

            @. x = [0.3, 0.4]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = [-5.0, 0.0]
            u = [0.3, 10.0]

        end

        fmin = Array{Function}(undef, m)

        for i=1:m

            fmin[i] = x -> ( 2.0 + 2.0 * i - ( exp( i * x[1] ) + exp( i * x[2] ) ) ) ^ 2.0

        end

    elseif nprob == 14

        # Brown and Dennis function

        @assert ( n == 4 ) && ( m ≥ n ) "The dimensions of the problem must satisfy n = 4 e m ≥ n."

        if usual_init_point

            @. x = [25.0, 5.0, - 5.0, - 1.0]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = [0.0, 0.0, -10.0, -10.0]
            u = [35.0, 10.0, 10.0, 0.0]

        end

        fmin = Array{Function}(undef, m)

        for i = 1:m

            ti = i / 5.0

            fmin[i] = x -> ( ( x[1] + ti * x[2] - exp( ti ) ) ^ 2.0 + ( x[3] + x[4] * sin( ti ) - cos( ti ) ) ^ 2.0 ) ^ 2.0

        end

    elseif nprob == 15

        # Chebyquad function

        @assert ( n ≥ 1 ) && ( m ≥ n ) "The dimensions of the problem must satisfy n ≥ 1 e m ≥ n."

        if usual_init_point

            c = n + 1.0

            for i=1:n

                x[i] = i / c

            end

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = zeros(Float64, n)
            l[n] = 0.75
            u = ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        function f15a(x, i)

            y = 2.0 * x - 1.0

            if y > 1.0

                return cosh( i * acosh( y ) )

            elseif y < - 1.0

                return ( - 1.0 ) ^ i * cosh( i * acosh( - y ) )

            else

                return cos( i * acos( y ) )

            end

        end

        function f15b(x, i)

            n = length(x)

            s = 0.0

            for j=1:n

                a = f15a( x[j], i )
                s += a

            end

            s *= ( 1.0 / n )

            if iseven( i )

                s += 1.0 / ( i ^ 2.0 - 1.0 )

            end

            return s

        end

        for i=1:m

            fmin[i] = x -> ( f15b(x, i) ) ^ 2.0

        end

    elseif nprob == 16

        # Brown almost-linear function

        @assert ( n ≥ 1 ) && ( m == n ) "The dimensions of the problem must satisfy n ≥ 1 e m = n."

        if usual_init_point

            @. x = 0.5

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - ones(Float64, n)
            u = 0.75 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        for i = 1:(n-1)

            fmin[i] = x -> ( x[i] + sum(x) - ( n + 1.0 ) ) ^ 2.0

        end

        fmin[n] = x -> ( prod(x) - 1.0 ) ^ 2.0

    elseif nprob == 17

        # Osborne 1 function

        @assert ( n == 5 ) && ( m == 33 ) "The dimensions of the problem must satisfy n = 5 e m = 33."

        if usual_init_point

            @. x = [0.5, 1.5, -1.0, 0.01, 0.02]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - ones(Float64, n)
            u = 2.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)
        y = [
                0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.850, 0.818, 0.784, 0.751,
                0.718, 0.685, 0.658, 0.628, 0.603, 0.580, 0.558, 0.538, 0.522, 0.506, 0.490,
                0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.420, 0.414, 0.411, 0.406
                ]

        for i = 1:m

            ti = 10.0 * (i - 1.0)
            fmin[i] = x -> ( y[i] - ( x[1] + x[2] * exp( - ti * x[4] ) + x[3] * exp( - ti * x[5] ) ) ) ^ 2.0

        end

    elseif nprob == 18

        # Osborne 2 function

        @assert ( n == 11 ) && ( m == 65 ) "The dimensions of the problem must satisfy n = 11 e m = 65."

        if usual_init_point

            @. x = [1.3, 0.65, 0.65, 0.7, 0.6, 3.0, 5.0, 7.0, 2.0, 4.5, 5.5]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = zeros(Float64, n)
            l[n] = 4.0
            u = 50.0 * ones(Float64, n)
            u[1] = 14.0
            u[2] = 5.0
            u[3] = 5.0
            u[4] = 5.0

        end

        fmin = Array{Function}(undef, m)
        y = [
                1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725, 0.746, 0.679, 0.608,
                0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724, 0.649, 0.649, 0.694, 0.644, 0.624, 0.661,
                0.612, 0.558, 0.533, 0.495, 0.500, 0.423, 0.395, 0.375, 0.372, 0.391, 0.396, 0.405, 0.428,
                0.429, 0.523, 0.562, 0.607, 0.653, 0.672, 0.708, 0.633, 0.668, 0.645, 0.632, 0.591, 0.559,
                0.597, 0.625, 0.739, 0.710, 0.729, 0.720, 0.636, 0.581, 0.428, 0.292, 0.162, 0.098, 0.054
                ]

        for i=1:m

            ti = (i - 1.0) / 10.0
            fmin[i] = x -> ( y[i] - ( x[1] * exp( - ti * x[5] ) + x[2] * exp( - ( ti - x[9] ) ^ 2.0 * x[6] ) + x[3] * exp( - ( ti - x[10] ) ^ 2.0 * x[7] ) + x[4] * exp( - ( ti - x[11] ) ^ 2.0 * x[8] ) ) ) ^ 2.0

        end

    elseif nprob == 19

        # BDQRTIC function

        @assert ( n ≥ 5 ) && ( m == 2 * ( n - 4 ) ) "The dimensions of the problem must satisfy n ≥ 5 e m = 2 * ( n - 4 )."

        if usual_init_point

            @. x = 1.0

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 5.0 * ones(Float64, n)
            l[1] = 0.8
            u = 10.0 * ones(Float64, n)

        end

        fmin = Array{Function}(undef, m)

        for i=1:(n - 4)

            fmin[i] = x -> ( - 4.0 * x[i] + 3.0 ) ^ 2.0
            fmin[n - 4 + i] = x -> ( x[i] ^ 2.0 + 2.0 * x[i + 1] ^ 2.0 + 3.0 * x[i + 2] ^ 2.0 + 4.0 * x[i + 3] ^ 2.0 + 5.0 * x[n] ^ 2.0 ) ^ 2.0

        end

    elseif nprob == 20

        # Cube function

        @assert ( n ≥ 2 ) && ( m == n ) "The dimensions of the problem must satisfy n ≥ 2 e m = n."

        if usual_init_point

            @. x = 0.5

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 5.0 * ones(Float64, n)
            u = 10.0 * ones(Float64, n)
            u[1] = 0.8

        end

        fmin = Array{Function}(undef, m)

        fmin[1] = x -> ( x[1] - 1.0 ) ^ 2.0

        for i=2:m

            fmin[i] = x -> ( 10.0 * ( x[i] - x[i - 1] ^ 3.0 ) ) ^ 2.0

        end

    elseif nprob == 21

        # Mancino function

        @assert ( n ≥ 2 ) && ( m == n ) "The dimensions of the problem must satisfy n ≥ 2 e m = n."

        if usual_init_point

            for i = 1:n

                s = 0.0

                for j = 1:n

                    ti = i / j

                    s += sqrt( ti ) * ( ( sin( log( sqrt( ti ) ) ) ) ^ 5.0 + ( cos( log( sqrt( ti ) ) ) ) ^ 5.0  )

                end

                x[i] = - 8.710996 - 4.0 * ( ( i - 50.0 ) ^ 3.0 + s )

            end

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = 2.0e5 * ones(Float64, n)
            l[n-1] = 2.5e5
            u = 1.0e7 * ones(Float64, n)
            u[n] = 3.0e5

        end

        fmin = Array{Function}(undef, m)

        function f21(x, i)

            s = 0.0

            for j = 1:length(x)

                v = sqrt( x[i] ^ 2.0 + i / j )
                s += v * ( ( sin( log( v ) ) ) ^ 5.0 + ( cos( log( v ) ) ) ^ 5.0 )

            end

            return s

        end

        for i = 1:m

            fmin[i] = x -> ( 1400.0 * x[i] + ( i - 50.0 ) ^ 3.0 + f21(x, i) ) ^ 2.0

        end

    elseif nprob == 22

        # HEART8LS function

        @assert ( n == 8 ) && ( m == 8 ) "The dimensions of the problem must satisfy n = m = 8."

        if usual_init_point

            @. x = [ - 0.3, - 0.39, 0.3, - 0.344, - 1.2, 2.69, 1.59, - 1.5 ]

        end

        if unconstrained

            l = c_lower * ones(Float64, n)
            u = c_upper * ones(Float64, n)

        else

            l = - 30.0 * ones(Float64, n)
            u = 30.0 * ones(Float64, n)
            u[3] = 3.5
            u[n] = -1.0

        end

        fmin = Array{Function}(undef, m)
        fmin[1] = x -> ( x[1] + x[2] + 0.69 ) ^ 2.0
        fmin[2] = x -> ( x[3] + x[4] + 0.044 ) ^ 2.0
        fmin[3] = x -> ( x[5] * x[1] + x[6] * x[2] - x[7] * x[3] - x[8] * x[4] + 1.57 ) ^ 2.0
        fmin[4] = x -> ( x[7] * x[1] + x[8] * x[2] + x[5] * x[3] + x[6] * x[4] + 1.31 ) ^ 2.0
        fmin[5] = x -> ( x[1] * ( x[5] ^ 2.0 - x[7] ^ 2.0 ) - 2.0 * x[3] * x[5] * x[7] + x[2] * ( x[6] ^ 2.0 - x[8] ^ 2.0 ) - 2.0 * x[4] * x[6] * x[8] + 2.65 ) ^ 2.0
        fmin[6] = x -> ( x[3] * ( x[5] ^ 2.0 - x[7] ^ 2.0 ) + 2.0 * x[1] * x[5] * x[7] + x[4] * ( x[6] ^ 2.0 - x[8] ^ 2.0 ) + 2.0 * x[2] * x[6] * x[8] - 2.0 ) ^ 2.0
        fmin[7] = x -> ( x[1] * x[5] * ( x[5] ^ 2.0 - 3.0 * x[7] ^ 2.0 ) + x[3] * x[7] * ( x[7] ^ 2.0 - 3.0 * x[5] ^ 2.0 ) + x[2] * x[6] * ( x[6] ^ 2.0 - 3.0 * x[8] ^ 2.0 ) + x[4] * x[8] * ( x[8] ^ 2.0 - 3.0 * x[6] ^ 2.0 ) + 12.6 ) ^ 2.0
        fmin[8] = x -> ( x[3] * x[5] * ( x[5] ^ 2.0 - 3.0 * x[7] ^ 2.0 ) - x[1] * x[7] * ( x[7] ^ 2.0 - 3.0 * x[5] ^ 2.0 ) + x[4] * x[6] * ( x[6] ^ 2.0 - 3.0 * x[8] ^ 2.0 ) - x[2] * x[8] * ( x[8] ^ 2.0 - 3.0 * x[6] ^ 2.0 ) - 9.48 ) ^ 2.0

    end

    if rsp

        @. x *= 10.0

    end

    return x, l, u, fmin

end

function problem_generator_hs(
                                nprob::Int64
                                )

    @assert 1 ≤ nprob ≤ 100 "The number of the problem must satisfy 1 ≤ nprob ≤ 100."

    # Defines the objective functions based on Hock-Schittkowski test-set.

    function hs1(x)

        return 100.0 * ( x[2] - x[1] ^ 2.0 ) ^ 2.0 + ( 1.0 - x[1] ) ^ 2.0

    end

    function hs3(x)

        return x[2] + 1.0e-5 * ( x[2] - x[1] ) ^ 2.0

    end

    function hs4(x)

        return ( x[1]  + 1.0 ) ^ 3.0 / 3.0 + x[2]

    end

    function hs5(x)

        return sin( x[1] + x[2] ) + ( x[1] - x[2] ) ^ 2.0 - 1.5 * x[1] + 2.5 * x[2] + 1.0

    end

    function hs25(x)

        s = 0.0

        for i = 1:99

            ui = 25.0 + ( - 50.0 * log( 0.01 * i ) ) ^ ( 2.0 / 3.0)

            s += ( - 0.01 * i + exp( - ( 1.0 / x[1] ) * ( ui - x[2] ) ^ ( x[3] ) ) ) ^ 2.0

        end

        return s

    end

    function hs38(x)

        return 100.0 * ( x[2] - x[1] ^ 2.0 ) ^ 2.0 + ( 1.0 - x[1] ) ^ 2.0 + 90.0 * ( x[4] - x[3] ^ 2.0 ) ^ 2.0 + ( 1.0 - x[3] ) ^ 2.0 + 10.1 * ( ( x[2] - 1.0 ) ^ 2.0 + ( x[4] - 1.0 ) ^ 2.0 ) + 19.8 * ( x[2] - 1.0 ) * ( x[4] - 1.0 )

    end

    function hs45(x)

        return 2.0 - x[1] * x[2] * x[3] * x[4] * x[5] / 120.0

    end

    function hs110(x)

        s1 = 0.0
        s2 = 1.0

        for i = 1:10

            s1 += ( log( x[i] - 2.0 ) ) ^ 2.0 + ( log( 10.0 - x[i] ) ) ^ 2.0
            s2 *= x[i]

        end

        return s1 - s2 ^ ( 0.2 )

    end

    # Initializes some useful variables and constants
    c_lower = - 1.0e20
    c_upper = 1.0e20

    if nprob == 1

        n = 2
        x = [ -2.0, 1.0 ]
        l = [ c_lower, 0.0 ]
        u = c_upper * ones(Float64, n)
        fmin = [ hs1, hs3 ]

    elseif nprob == 2

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = c_upper * ones(Float64, n)
        fmin = [ hs1, hs4 ]

    elseif nprob == 3

        n = 2
        x = zeros(Float64, n)
        l = [ - 1.5, - 1.5 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs1, hs5 ]

    elseif nprob == 4

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 0.1, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs1, hs25 ]

    elseif nprob == 5

        n = 4
        x = [ - 3.0, - 1.0, - 3.0, - 1.0 ]
        l = [ - 10.0, - 1.5, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, n)
        fmin = [ hs1, hs38 ]

    elseif nprob == 6

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs45 ]

    elseif nprob == 7

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs110 ]

    elseif nprob == 8

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = c_upper * ones(Float64, n)
        fmin = [ hs3, hs4 ]

    elseif nprob == 9

        n = 2
        x = zeros(Float64, n)
        l = [ - 1.5, 0.0 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs3, hs5 ]

    elseif nprob == 10

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 0.1, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs3, hs25 ]

    elseif  nprob == 11

        n = 4
        x = [ - 3.0, - 1.0, - 3.0, - 1.0 ]
        l = [ - 10.0, 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, 4)
        fmin = [ hs3, hs38 ]

    elseif nprob == 12

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs45 ]

    elseif nprob == 13

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs3, hs110 ]

    elseif nprob == 14

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs4, hs5 ]

    elseif nprob == 15

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs4, hs25 ]

    elseif nprob == 16

        n = 4
        x = [ 3.0, 1.0, - 3.0, - 1.0 ]
        l = [ 1.0 , 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, n)
        fmin = [ hs4, hs38 ]

    elseif nprob == 17

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs45 ]

    elseif nprob == 18

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs4, hs110 ]

    elseif nprob == 19

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ -1.5, -3.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs5, hs38 ]

    elseif nprob == 20

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs5, hs45 ]

    elseif nprob == 21

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs25, hs38 ]

    elseif nprob == 22

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs25, hs45 ]

    elseif nprob == 23

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs38, hs45 ]

    elseif nprob == 24

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs38, hs110 ]

    elseif nprob == 25

        n = 2
        x = [10.0, 1.0]
        l = [1.0, 0.0]
        u = c_upper * ones(Float64, n)
        fmin = [ hs1, hs3, hs4 ]

    elseif nprob == 26

        n = 2
        x = [10.0, 1.0]
        l = [- 1.5, 0.0]
        u = [4.0, 3.0]
        fmin = [ hs1, hs3, hs5 ]

    elseif nprob == 27

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 0.1, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs1, hs3, hs25 ]

    elseif  nprob == 28

        n = 4
        x = [ - 3.0, - 1.0, - 3.0, - 1.0 ]
        l = [ - 10.0, 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, 4)
        fmin = [ hs1, hs3, hs38 ]

    elseif nprob == 29

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs3, hs45 ]

    elseif nprob == 30

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs3, hs110 ]

    elseif nprob == 31

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs1, hs4, hs5 ]

    elseif nprob == 32

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs1, hs4, hs25 ]

    elseif nprob == 33

        n = 4
        x = [ 3.0, 1.0, - 3.0, - 1.0 ]
        l = [ 1.0 , 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, n)
        fmin = [ hs1, hs4, hs38 ]

    elseif nprob == 34

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs4, hs45 ]

    elseif nprob == 35

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs4, hs110 ]

    elseif nprob == 36

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ -1.5, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs1, hs5, hs38 ]

    elseif nprob == 37

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs5, hs45 ]

    elseif nprob == 38

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs1, hs25, hs38 ]

    elseif nprob == 39

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs25, hs45 ]

    elseif nprob == 40

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs38, hs45 ]

    elseif nprob == 41

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs38, hs110 ]

    elseif nprob == 42

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs3, hs4, hs5 ]

    elseif nprob == 43

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs3, hs4, hs25 ]

    elseif nprob == 44

        n = 4
        x = [ 3.0, 1.0, - 3.0, - 1.0 ]
        l = [ 1.0 , 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, n)
        fmin = [ hs3, hs4, hs38 ]

    elseif nprob == 45

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs4, hs45 ]

    elseif nprob == 46

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs3, hs4, hs110 ]

    elseif nprob == 47

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ -1.5, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs3, hs5, hs38 ]

    elseif nprob == 48

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs5, hs45 ]

    elseif nprob == 49

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs3, hs25, hs38 ]

    elseif nprob == 50

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs25, hs45 ]

    elseif nprob == 51

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs38, hs45 ]

    elseif nprob == 52

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs3, hs38, hs110 ]

    elseif nprob == 53

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 4.0, 3.0, 5.0 ]
        fmin = [ hs4, hs5, hs25 ]

    elseif nprob == 54

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs4, hs5, hs38 ]

    elseif nprob == 55

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs5, hs45 ]

    elseif nprob == 56

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs4, hs25, hs38 ]

    elseif nprob == 57

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs25, hs45 ]

    elseif nprob == 58

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs38, hs45 ]

    elseif nprob == 59

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs4, hs38, hs110 ]

    elseif nprob == 60

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 4.0, 3.0, 5.0, 10.0 ]
        fmin = [ hs5, hs25, hs38 ]

    elseif nprob == 61

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs5, hs25, hs45 ]

    elseif nprob == 62

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs5, hs38, hs45 ]

    elseif nprob == 63

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs25, hs38, hs45 ]

    elseif nprob == 64

        n = 2
        x = [ 1.125, 0.125 ]
        l = [ 1.0, 0.0 ]
        u = [ 4.0, 3.0 ]
        fmin = [ hs1, hs3, hs4, hs5 ]

    elseif nprob == 65

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 100.0, 25.6, 5.0 ]
        fmin = [ hs1, hs3, hs4, hs25 ]

    elseif nprob == 66

        n = 4
        x = [ 3.0, 1.0, - 3.0, - 1.0 ]
        l = [ 1.0 , 0.0, - 10.0, - 10.0 ]
        u = 10.0 * ones(Float64, n)
        fmin = [ hs1, hs3, hs4, hs38 ]

    elseif nprob == 67

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs3, hs4, hs45 ]

    elseif nprob == 68

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs3, hs4, hs110 ]

    elseif nprob == 69

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 0.1, 0.0, 0.0 ]
        u = [ 4.0, 3.0, 5.0 ]
        fmin = [ hs1, hs3, hs5, hs25 ]

    elseif nprob == 70

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ -1.5, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs1, hs3, hs5, hs38 ]

    elseif nprob == 71

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs3, hs5, hs45 ]

    elseif nprob == 72

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs1, hs3, hs25, hs38 ]

    elseif nprob == 73

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs3, hs25, hs45 ]

    elseif nprob == 74

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs3, hs38, hs45 ]

    elseif nprob == 75

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs3, hs38, hs110 ]

    elseif nprob == 76

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 4.0, 3.0, 5.0 ]
        fmin = [ hs1, hs4, hs5, hs25 ]

    elseif nprob == 77

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs1, hs4, hs5, hs38 ]

    elseif nprob == 78

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs4, hs5, hs45 ]

    elseif nprob == 79

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs1, hs4, hs25, hs38 ]

    elseif nprob == 80

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs4, hs25, hs45 ]

    elseif nprob == 81

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs4, hs38, hs45 ]

    elseif nprob == 82

        n = 10
        x = 9.0 * ones(Float64, n)
        l = 2.001 * ones(Float64, n)
        u = 9.999 * ones(Float64, n)
        fmin = [ hs1, hs4, hs38, hs110 ]

    elseif nprob == 83

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 4.0, 3.0, 5.0, 10.0 ]
        fmin = [ hs1, hs5, hs25, hs38 ]

    elseif nprob == 84

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs5, hs25, hs45 ]

    elseif nprob == 85

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs5, hs38, hs45 ]

    elseif nprob == 86

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs1, hs25, hs38, hs45 ]

    elseif nprob == 87

        n = 3
        x = [ 100.0, 12.5, 3.0]
        l = [ 1.0, 0.0, 0.0 ]
        u = [ 4.0, 3.0, 5.0 ]
        fmin = [ hs3, hs4, hs5, hs25 ]

    elseif nprob == 88

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, -10.0, -10.0 ]
        u = [ 4.0, 3.0, 10.0, 10.0 ]
        fmin = [ hs3, hs4, hs5, hs38 ]

    elseif nprob == 89

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs4, hs5, hs45 ]

    elseif nprob == 90

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, 0.0, -10.0 ]
        u = [ 10.0, 10.0, 5.0, 10.0 ]
        fmin = [ hs3, hs4, hs25, hs38 ]

    elseif nprob == 91

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs4, hs25, hs45 ]

    elseif nprob == 92

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 0.1, 0.0, 0.0, -10.0 ]
        u = [ 4.0, 3.0, 5.0, 10.0 ]
        fmin = [ hs3, hs5, hs25, hs38 ]

    elseif nprob == 93

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs5, hs25, hs45 ]

    elseif nprob == 94

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs5, hs38, hs45 ]

    elseif nprob == 95

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs3, hs25, hs38, hs45 ]

    elseif nprob == 96

        n = 4
        x = [ -3.0, -1.0, -3.0, -1.0 ]
        l = [ 1.0, 0.0, 0.0, -10.0 ]
        u = [ 4.0, 3.0, 5.0, 10.0 ]
        fmin = [ hs4, hs5, hs25, hs38 ]

    elseif nprob == 97

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs5, hs25, hs45 ]

    elseif nprob == 98

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs5, hs38, hs45 ]

    elseif nprob == 99

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 1.0
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs4, hs25, hs38, hs45 ]

    elseif nprob == 100

        n = 5
        x = 2.0 * ones(Float64, n)
        l = zeros(Float64, n)
        l[1] = 0.1
        u = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
        fmin = [ hs5, hs25, hs38, hs45 ]

    end

    return x, l, u, fmin

end


text:
module SymmetricSquares

export e, r90, r180, r270, tr, trr, fx, fy

include("symmetry.jl")

end

ext: jl


text:
@testset "Square root Augmented Unscented Kalman filter" begin
    @testset "time update" begin
        x = randn(6)
        A = randn(6,6)
        P = A'A
        P_chol = cholesky(P)
        B = randn(6,6)
        Q = B'B
        Q_chol = cholesky(Q)
        F = randn(6,6)
        f(x) = F * x
        f(x, noise) = F * x .+ noise

        tu = time_update(x, P, F, Q)
        tu_chol = time_update(x, P_chol, f, Augment(Q_chol))
        @test get_covariance(tu_chol) ≈ get_covariance(tu)
        @test get_state(tu_chol) ≈ get_state(tu)

        f!(y, x) = mul!(y, F, x)
        f!(y, x, noise) = y.= @~ F * x .+ noise
        tu_inter = SRAUKFTUIntermediate(6)
        tu_aug_inplace = time_update!(tu_inter, x, P_chol, f!, Augment(Q_chol))
        @test get_covariance(tu_aug_inplace) ≈ get_covariance(tu)
        @test get_state(tu_aug_inplace) ≈ get_state(tu)
    end

    @testset "measurement update" begin
        x = randn(6)
        A = randn(6,6)
        P = A'A
        P_chol = cholesky(P)
        B = randn(3,3)
        R = B'B
        y = randn(3)
        R_chol = cholesky(R)
        H = randn(3,6)
        h(x) = H * x
        h(x, noise) = H * x .+ noise

        mu = measurement_update(x, P, y, H, R)
        mu_chol = @inferred measurement_update(x, P_chol, y, h, Augment(R_chol))
        @test get_covariance(mu) ≈ get_covariance(mu_chol)
        @test get_state(mu) ≈ get_state(mu_chol)

        h!(y, x) = mul!(y, H, x)
        h!(y, x, noise) = y .= @~ H * x .+ noise
        mu_inter = @inferred SRAUKFMUIntermediate(6, 3)
        mu_ukf_inplace = measurement_update!(mu_inter, x, P_chol, y, h!, Augment(R_chol))
        @test get_covariance(mu_ukf_inplace) ≈ get_covariance(mu)
        @test get_state(mu_ukf_inplace) ≈ get_state(mu)
    end
end


text:
# Plots an random example of a prior and the resulting information distribution
#
# This code uses the planner dataset to produce more relevant results

using SubmodularMaximization
using JLD2
using Statistics
using Profile
using ProfileView

println("Loading rollout library")
library_file = "./data/rollout_library.jld2"
@time @load library_file data num_robots trials
rollouts = data

sparse = true

num_information_samples = 10
num_trials = 100

# Overwrite the number of robots
num_robots = 10

if sparse
  global copy_filter(x) = SparseFilter(x; threshold=1e-3)
else
  global copy_filter(x) = Filter(x)
end

# Full problem instances for each trial
trials_data = vcat(map(trials) do trial
                     rollout = rollouts[num_robots, trial]
                     configs = rollout.configs

                     [(grid = configs.grid,
                       sensor = configs.sensor,
                       target_filters = map(copy_filter, sample.target_filters),
                       trajectories = sample.trajectories,
                      )
                      for sample in rollout.samples]
                   end...)

num_filters = length(first(trials_data).target_filters)
samples_per_evaluation = num_filters * num_information_samples

# Create some fake observations and update the filter
function run_test(num_trials; print=true)
  ifprint(xs...) = if print println(xs...) end

  time = @elapsed for trial = trials_data[1:num_trials]
    mean(trial.target_filters) do filter

      finite_horizon_information(trial.grid,
                                 filter,
                                 trial.sensor,
                                 trial.trajectories,
                                 num_samples = num_information_samples
                                ).reward
    end
  end

  ifprint("Test duration: ", time, " seconds")
  ifprint("  ", time / num_trials, " seconds per information estimate")
  ifprint("  ", time / num_trials / samples_per_evaluation,
          " seconds per information estimate per sample")
  ifprint("Note: total number of information samples is ",
          samples_per_evaluation)
end

run_test(1)
@time run_test(num_trials)

## Run the profiler
@profview run_test(1, print=false)

print("Waiting for input...")
readline()

#ProfileView.closeall()
@profview run_test(num_trials, print=false)

println()
Profile.print(mincount=100, format=:flat, sortedby=:count)

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
import TrajectoryOptimization: primals, duals, packZ, update!
using BlockArrays

# Create solver
opts = ProjectedNewtonSolverOptions{Float64}()
solver = SequentialNewtonSolver(prob,opts)
NN = N*n + (N-1)*m

# Test update functions
TO.dynamics_constraints!(prob, solver)
TO.dynamics_jacobian!(prob, solver)
TO.update_constraints!(prob, solver)
TO.constraint_jacobian!(prob, solver)
TO.active_set!(prob, solver)
TO.cost_expansion!(prob, solver)
TO.invert_hessian!(prob, solver)

# Init solvers
solver0 = ProjectedNewtonSolver(prob,opts)
TO.update!(prob, solver0)
solver = SequentialNewtonSolver(prob, opts)
TO.update!(prob, solver)

# Projection
Hinv = inv(Diagonal(solver0.H))
Y,y0 = TO.active_constraints(prob, solver0)
S0 = Y*Hinv*Y'
δz0 = -Hinv*Y'*(S0\y0)

y = TO.active_constraints(solver)
vcat(y...) ≈ y0
TO.calc_factors!(solver)

L = Array(TO.buildL(solver))
L - cholesky(Symmetric(Array(S0))).L

a = solver.active_set
S = TO.buildS(solver,C,D)
S = zero(S0)
S = zero(Array(S))
TO.buildS!(S,solver,C,D)
@btime TO.buildS!($S,$solver,$C,$D)
@btime TO.buildS($solver,$C,$D)
@btime $Y*$Hinv*$Y'

Hinv = Diagonal(Array(Hinv))
nnz(Y)/length(Y)
Y = Array(Y)

S ≈ S0
S = TO.buildShurCompliment(prob, solver)
@btime TO.buildShurCompliment($prob, $solver)

y_part = TO.dual_partition(solver)
Pa = sum(y_part)
S = BlockArray(zeros(Pa,Pa),y_part,y_part)
S[Block(2,2)]
view(S,Block(2,2)).indices

S_part = TO.Sinds(solver)
S = zeros(Pa,Pa)
using Profile
Profile.init(n=10000000,delay=1e-9)
@profiler [TO.buildS!(S,solver,A,B,C,D,S_part) for k = 1:100]
@code_warntype TO.buildS!(S,solver,A,B,C,D,S_part)
@which TO.buildS!(S,solver,A,B,C,D,S_part)
Symmetric(S) ≈ S0
A = [solver.∇F[k].xx for k = 2:N]
B = [solver.∇F[k].xu for k = 2:N]

@btime TO.buildS!($S,$solver,$A,$B,$C,$D,$S_part)
@btime $Y*$Hinv*$Y'

@btime cholesky(Symmetric($S0))

S0_ = BlockArray(S0,y_part,y_part)

# Multiplier projection
λ0 = duals(solver0.V)[solver0.a.duals]
res0 = solver0.g + Y'λ0

λa = TO.active_duals(solver.V, solver.active_set)
δλ = TO._mult_projection(solver, λa)
@test vcat(δλ...) ≈ -(Y*Y')\(Y*res0)

# KKT solve
S0 = Symmetric(Y*Hinv*Y')
r = solver0.g + Y'λ0
δλ0 = cholesky(S0)\(y0 - Y*Hinv*r)
δz0 = -Hinv*(r + Y'δλ0)
δV1,δλ1 = TO.solveKKT_Shur(prob,solver0,Hinv)

δx, δu, δλ = TO._solveKKT(solver, λa)
@test packZ(δx, δu) ≈ δz0
@test packZ(δx, δu) ≈ δV1[1:NN]
@test vcat(δλ...) ≈ δλ0
@test vcat(δλ...) ≈ δλ1

@btime TO.solveKKT_Shur($prob,$solver0,$Hinv)

@btime TO.solveKKT($solver, $solver.V)

# Update methods
δV = zero(solver0.V.V)
δV[solver0.parts.primals] .= δz0
δV[solver0.parts.duals[solver0.a.duals]] .= δλ0
V_0 = solver0.V + δV

V_0 = solver0.V + δV
V_ = copy(solver.V)
TO._update_primals!(V_, δx, δu)
TO._update_duals!(V_, δλ, solver.active_set)
@test V_.X ≈ V_0.X
@test V_.U ≈ V_0.U
@test vcat(V_.λ...) ≈ V_0.Y

V2_ = solver.V + (δx, δu)
@test V2_.X == V_.X
@test V2_.U == V_.U
@test V2_.λ == solver.V.λ

V2_ = solver.V + (δλ, solver.active_set)
@test V2_.λ == V_.λ

V2_ = solver.V + ((δx, δu, δλ), solver.active_set)
@test V2_.X == V_.X
@test V2_.U == V_.U
@test V2_.λ == V_.λ

# Newton Step from beginning
opts = ProjectedNewtonSolverOptions{Float64}()
opts.verbose = true
opts.feasibility_tolerance = 1e-10
opts.active_set_tolerance = 1e-4
solver0 = ProjectedNewtonSolver(prob, opts)
update!(prob, solver0)
solver = SequentialNewtonSolver(prob, opts)
update!(prob, solver)
V = solver.V

# Projection
TO.projection!(prob, solver0)
TO.projection!(prob, solver, V)

# Multiplier projection
norm(TO.residual(prob,solver0))
r0, δλ0 = TO.multiplier_projection!(prob, solver0)
TO.res2(solver, V)
r, δλ = TO.multiplier_projection!(solver, V)
@test r ≈ r0
@test vcat(δλ...) ≈ δλ0

# KKT Solve
J0 = cost(prob, solver0.V)
viol0 = max_violation(solver0)
Y,y0 = TO.active_constraints(prob, solver0)
λ0 = duals(solver0.V)[solver0.a.duals]
δV0,δλ0 = TO.solveKKT_Shur(prob, solver0, Hinv)

J = cost(prob, V)
viol = max_violation(solver)
λa = TO.active_duals(V, solver.active_set)
δx, δu, δλ = TO.solveKKT(solver, V)

@test J0 ≈ J
@test viol0 ≈ viol
@test packZ(δx, δu) ≈ δV0[1:NN]
@test norm(vcat(δλ...) - δλ0) < 1e-10

res_init = TO.res2(solver,V)
@test norm(TO.residual(prob,solver0)) ≈ res_init

# Line Search
TO.update!(prob, solver0)
TO.update!(prob, solver)
α = 0.9
V_0 = solver0.V + α*δV0
V_ = V + (α.*(δx,δu,δλ),solver.active_set)

@test V_0.X ≈ V_.X
@test V_0.U ≈ V_.U
@test V_0.Y ≈ vcat(V_.λ...)

TO.dynamics_constraints!(prob, solver0, V_0)
TO.update_constraints!(prob, solver0, V_0)
TO.dynamics_jacobian!(prob, solver0, V_0)
TO.constraint_jacobian!(prob, solver0, V_0)
TO.active_set!(prob, solver0)
Y0,y0 = TO.active_constraints(prob,solver0)
viol0 = norm(y,Inf)
S0 = Y0*Hinv*Y0'
δZ = -Hinv*Y0'*((Y0*Hinv*Y0')\y0)
δλ0 = (Y0*Hinv*Y0')\y0

TO.dynamics_constraints!(prob, solver, V_)
TO.update_constraints!(prob, solver, V_)
TO.dynamics_jacobian!(prob, solver, V_)
TO.constraint_jacobian!(prob, solver, V_)
TO.active_set!(prob, solver)
TO.cost_expansion!(prob, solver, V_)
y = TO.active_constraints(solver)
maximum(norm.(y,Inf))
x,u = TO._projection(solver)

@btime TO.calc_factors!($solver)
@btime cholesky(Symmetric(S0))
@btime
TO.calc_factors!(solver)
δλ = TO.solve_cholesky(solver,y)
vcat(δλ...) ≈ δλ0
vcat(y...) ≈ y0
packZ(x,u) ≈ δZ

Y = TO.buildY(solver)
Y ≈ Y0
L = TO.buildL(solver)

L0 = cholesky(Symmetric(Array(S0))).L
L0 = PseudoBlockArray(Array(L0),y_part,y_part)

D = L0 .≈ Array(L)
D = PseudoBlockArray(Array(D),y_part,y_part)
D[Block(6,5)]

norm(Array(L0) - Array(L),Inf)
L0[Block(2N,2N)] ≈ L[Block(6,5)]

y_part = TO.dual_partition(solver)
[all(D[Block(i,j)]) for i=1:2N,j=1:2N]

L0 = PseudoBlockArray(L,y_part,y_part)
L0[Block(6,5)] - L[Block(6,5)]

all(D[Block(6,5)])
[all(D[Block(k,k)]) for k = 1:2N-2]
Array{Int}(L*L' .≈ (Y*Hinv*Y'))

TO.projection!(prob, solver0, V_0)
TO.projection!(prob, solver, V_)
res0, = TO.multiplier_projection!(prob, solver0, V_0)
res, = TO.multiplier_projection!(solver, V_)
@test res0 ≈ res
J0 = cost(prob, V_0)
J = cost(prob, V_)
@test J0 ≈ J

TO.line_search(prob, solver0, δV0)
TO.line_search(prob, solver, δx, δu, δλ)

# Test Whole Step
opts = ProjectedNewtonSolverOptions{Float64}(feasibility_tolerance=1e-10)
opts.verbose = true
opts.active_set_tolerance = 1e-4

solver0 = ProjectedNewtonSolver(prob,opts)
begin
    copyto!.(solver0.V.X, prob.X)
    copyto!.(solver0.V.U, prob.U)
    TO.newton_step!(prob, solver0)
end

solver = SequentialNewtonSolver(prob, opts)
begin
    copyto!(solver.V.X, prob.X)
    copyto!(solver.V.U, prob.U)
    TO.newton_step!(prob, solver)
end

solver = SequentialNewtonSolver(prob, opts)
V_ = TO.newton_step!(prob, solver)
V_.X

J0 = cost(prob)
viol0 = max_violation(prob)
res = solve(prob,solver)
@test cost(res) < J0
@test max_violation(res) < viol0
@test max_violation(res) < 1e-10

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
s = Benchmarks.Samples()

@test isa(s.elapsed_times, Vector{Float64})
@test length(s.elapsed_times) == 0

@test isa(s.bytes_allocated, Vector{Int})
@test length(s.bytes_allocated) == 0

@test isa(s.gc_times, Vector{Float64})
@test length(s.gc_times) == 0

@test isa(s.num_allocations, Vector{Int})
@test length(s.num_allocations) == 0

e = Benchmarks.Environment()


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const unix = 1
const linux = 1

# begin enum ANONYMOUS_7
const ANONYMOUS_7 = UInt32
const TT_EMPTY = (UInt32)(0)
const TT_ERROR = (UInt32)(1)
const TT_LONG = (UInt32)(2)
const TT_DOUBLE = (UInt32)(3)
const TT_STRING = (UInt32)(4)
# end enum ANONYMOUS_7

# begin enum ANONYMOUS_8
const ANONYMOUS_8 = Cint
const VR_OK = (Int32)(0)
const VR_OUTOFMEMORY = (Int32)(-1)
const VR_BADVARTYPE = (Int32)(-2)
const VR_INVALIDARG = (Int32)(-3)
const VR_INVALIDROW = (Int32)(-4)
const VR_INVALIDCOL = (Int32)(-5)
# end enum ANONYMOUS_8

# begin enum ANONYMOUS_6
const ANONYMOUS_6 = Cint
const IPQ_OK = (Int32)(0)
const IPQ_OUTOFMEMORY = (Int32)(-1)
const IPQ_BADVARTYPE = (Int32)(-2)
const IPQ_INVALIDARG = (Int32)(-3)
const IPQ_INVALIDROW = (Int32)(-4)
const IPQ_INVALIDCOL = (Int32)(-5)
const IPQ_BADINSTANCE = (Int32)(-6)
# end enum ANONYMOUS_6

# begin enum ANONYMOUS_4
const ANONYMOUS_4 = UInt32
const TT_EMPTY = (UInt32)(0)
const TT_ERROR = (UInt32)(1)
const TT_LONG = (UInt32)(2)
const TT_DOUBLE = (UInt32)(3)
const TT_STRING = (UInt32)(4)
# end enum ANONYMOUS_4

# begin enum ANONYMOUS_5
const ANONYMOUS_5 = Cint
const VR_OK = (Int32)(0)
const VR_OUTOFMEMORY = (Int32)(-1)
const VR_BADVARTYPE = (Int32)(-2)
const VR_INVALIDARG = (Int32)(-3)
const VR_INVALIDROW = (Int32)(-4)
const VR_INVALIDCOL = (Int32)(-5)
# end enum ANONYMOUS_5


text:
const mainmeta = [
        "ageMonths",
        "age0to3mo",
        "age3to6mo",
        "age6to12mo",
        "age12moplus",
        "mother_HHS_Education",
        "race",
        "cogScore",
        "has_segmentation",
        "ECHOTPCoded",
        "assessmentDate",
        "ed",
        "rce"
    ]

const brainmeta = [
            #  "CortexVol",
            #  "CorticalWhiteMatterVol",
            #  "SubCortGrayVol",
            #  "TotalGrayVol",
            #  "BrainSegVol-to-eTIV",
            #  "CerebralWhiteMatterVol",
            #  "lhCorticalWhiteMatterVol",
            #  "lhCerebralWhiteMatterVol",
            #  "lhCortexVol",
                "Left-Thalamus",
                "Left-Lateral-Ventricle",
                "Left-Cerebellum-White-Matter",
                "Left-Cerebellum-Cortex",
                "Left-Caudate",
                "Left-Putamen",
                "Left-Pallidum",
                "Left-Hippocampus",
                "Left-Amygdala",
                "Left-Accumbens-area",
                "Left-VentralDC",
                "Left-choroid-plexus",
            #  "rhCorticalWhiteMatterVol",
            #  "rhCerebralWhiteMatterVol",
            #  "rhCortexVol",
                "Right-Thalamus",
                "Right-Lateral-Ventricle",
                "Right-Cerebellum-White-Matter",
                "Right-Cerebellum-Cortex",
                "Right-Caudate",
                "Right-Putamen",
                "Right-Pallidum",
                "Right-Hippocampus",
                "Right-Amygdala",
                "Right-Accumbens-area",
                "Right-VentralDC",
                "Right-choroid-plexus",
                "Brain-Stem",
                "CSF"
    ]

"""
Generates some data objects useful in many scripts.
Requires files the following files,
generated in notebooks 1-3:

- "data/wrangled/omnisamples.csv"
- "data/wrangled/etohsamples.csv"
- "data/wrangled/timepoints.csv"
- "data/wrangled/metabolites.csv"
- "data/wrangled/species.csv"
"""
function startup(; dfs=[:omni, :etoh, :tps, :complete_brain, :metabolites, :species])
    omni = CSV.read(joinpath(@__DIR__, "..", "data/wrangled/omnisamples.csv"), DataFrame)
    etoh = CSV.read(joinpath(@__DIR__, "..", "data/wrangled/etohsamples.csv"), DataFrame)

    tps, complete_brain = _gentps()

    metabolites = _genmetabolites(etoh, tps)

   species = _genspecies(omni, tps)

   return (; omni, etoh, tps, complete_brain, metabolites, species)[dfs]
end

function _gentps()
    tps  = CSV.read(joinpath(@__DIR__, "..", "data/wrangled/timepoints.csv"), DataFrame)

    DataFrames.transform!(groupby(tps, :subject), :mother_HHS_Education => (r->coalesce(r...)) => :hhs)
    tps.ed = categorical(tps.hhs; levels=[-8 , 2:7...], ordered=true)
    tps.ed = recode(tps.ed,
        -8 => missing,
        2 => "Junior high school",
        3 => "Some high school",
        4 => "High school grad",
        5 => "Some college",
        6 => "College grad",
        7 => "Grad/professional school")

    DataFrames.transform!(groupby(tps, :subject), :race => (r->coalesce(r...)) => :race)
    tps.rce = categorical(tps.race; ordered=true)
    tps.rce = recode(tps.rce,
        "American Indian or Alaska Native"=> "Other",
        "Some other race"                 => "Other",
        "Asian Indian"                    => "Asian",
        "Black or African American"       => "Black",
        missing                           => "Unknown"
    )
    levels!(tps.rce, ["White","Black","Asian","Mixed","Other","Unknown"])

    tps."Left-Thalamus" = map(eachrow(tps)) do row
        (t, p) = (row."Left-Thalamus", row."Left-Thalamus-Proper")
        all(ismissing, (t,p)) && return missing
        return max(coalesce(t, 0), coalesce(p, 0))
    end

    tps."Right-Thalamus" = map(eachrow(tps)) do row
        (t, p) = (row."Right-Thalamus", row."Right-Thalamus-Proper")
        all(ismissing, (t,p)) && return missing
        return max(coalesce(t, 0), coalesce(p, 0))
    end
    brainmeta = Resonance.brainmeta
    for m in brainmeta
        tps[!, m] .= tps[!, m] ./ map(x-> ismissing(x) || x == 0 ? missing : x, tps."EstimatedTotalIntraCranialVol")
    end

    select!(tps, ["subject", "timepoint", mainmeta..., brainmeta...])
    rename!(tps, Dict(k=> replace(k, "-"=>"_") for k in brainmeta))
    brainmeta = map(i-> replace(brainmeta[i], "-"=>"_"), eachindex(brainmeta))

    complete_brain = completecases(tps[:, brainmeta])
    return tps, complete_brain
end

function _genmetabolites(etoh_samples, tps)
    metabolites = CSV.read(joinpath(@__DIR__, "..", "data/wrangled/metabolites.csv"), DataFrame)
    ms = [Resonance.Metabolite(row[:uid], row[:Metabolite], row[:MZ], row[:RT]) for row in eachrow(metabolites)]
    metabolites = CommunityProfile(Matrix(metabolites[!, 9:end]), ms, MicrobiomeSample.(names(metabolites)[9:end]))
    set!(metabolites, leftjoin(etoh_samples, tps[!, ["subject",  "timepoint", mainmeta...]], on=[:subject, :timepoint], makeunique=true))

    return metabolites
end

function _genspecies(omni_samples, tps)
    species = CSV.read(joinpath(@__DIR__, "..", "data/wrangled/species.csv"), DataFrame)
    species = CommunityProfile(Matrix(species[!, 2:end]), taxon.(species[!, 1]), MicrobiomeSample.(names(species)[2:end]))
    set!(species, leftjoin(omni_samples, tps[!, ["subject", "timepoint", mainmeta...]], on=[:subject, :timepoint], makeunique=true))
    species = species[:, map(!ismissing, get(species, :subject)) .& map(!ismissing, get(species, :timepoint))]

    return species
end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
using Pkg
using SafeTestsets
const LONGER_TESTS = false

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

@time begin

if GROUP == "All" || GROUP == "Core"
  @time @safetestset "FiniteDiff Standard Tests" begin include("finitedifftests.jl") end
  @time @safetestset "Color Differentiation Tests" begin include("coloring_tests.jl") end
  @time @safetestset "Out of Place Tests" begin include("out_of_place_tests.jl") end
end

if GROUP == "All" || GROUP == "Downstream"
  activate_downstream_env()
  @time @safetestset "ODEs" begin
    import OrdinaryDiffEq
    include(joinpath(dirname(pathof(OrdinaryDiffEq)), "..", "test/interface/sparsediff_tests.jl"))
  end
end

end


text:
module RocketTapOnSubscribeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: tap_on_subscribe()" begin

    println("Testing: operator tap_on_subscribe()")

    run_proxyshowcheck("TapOnSubscribe", tap_on_subscribe(() -> begin end))

    @testset begin
        values = Int[]
        source = from(1:5) |> tap_on_subscribe(() -> push!(values, -1))
        actor  = keep(Int)

        @test values == []

        subscribe!(source, actor)

        @test actor.values == [ 1, 2, 3, 4, 5 ]
        @test values       == [ -1 ]

        subscribe!(source, actor)

        @test actor.values == [ 1, 2, 3, 4, 5, 1, 2, 3, 4, 5 ]
        @test values       == [ -1, -1 ]
    end

    @testset begin
        values = Int[]
        source = from(1:5) |> tap_on_subscribe(() -> push!(values, -1), TapBeforeSubscription())

        @test values == []

        subscribe!(source, (e) -> push!(values, e))

        @test values == [ -1, 1, 2, 3, 4, 5 ]

        subscribe!(source, (e) -> push!(values, e))

        @test values == [ -1, 1, 2, 3, 4, 5, -1, 1, 2, 3, 4, 5 ]
    end

    @testset begin
        values = Int[]
        source = from(1:5) |> tap_on_subscribe(() -> push!(values, -1), TapAfterSubscription())

        @test values == []

        subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1, 2, 3, 4, 5, -1 ]

        subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1, 2, 3, 4, 5, -1, 1, 2, 3, 4, 5, -1 ]
    end

    @testset begin
        values = Int[]
        source = completed(Int) |> tap_on_subscribe(() -> push!(values, -1))
        actor  = keep(Int)

        @test values == []

        subscribe!(source, actor)

        @test actor.values == [ ]
        @test values       == [ -1 ]
    end

end

end


text:
function logsumexp(arr::AbstractArray{T}) where {T <: Real}
    max_arr = maximum(arr)
    max_arr + log(sum(exp.(arr .- max_arr)))
end

function logsumexp(x1::Real, x2::Real)
    max_arr = max(x1, x2)
    max_arr + log(exp(x1 - max_arr) + exp(x2 - max_arr))
end

export logsumexp

########
# MCMC #
########

function mh(model::Generator, proposal::Generator, proposal_args_rest::Tuple, trace)
    model_args = get_call_record(trace).args
    proposal_args_forward = (trace, proposal_args_rest...,)
    forward_trace = simulate(proposal, proposal_args_forward)
    forward_score = get_call_record(forward_trace).score
    constraints = get_assignment(forward_trace)
    (new_trace, weight, discard) = update(
        model, model_args, NoChange(), trace, constraints)
    proposal_args_backward = (new_trace, proposal_args_rest...,)
    backward_trace = assess(proposal, proposal_args_backward, discard)
    backward_score = get_call_record(backward_trace).score
    if log(rand()) < weight - forward_score + backward_score
        # accept
        return new_trace
    else
        # reject
        return trace
    end
end

function mh(model::Generator, selector::SelectionFunction, selector_args::Tuple, trace)
    (selection, _) = select(selector, selector_args, get_assignment(trace))
    model_args = get_call_record(trace).args
    (new_trace, weight) = regenerate(model, model_args, NoChange(), trace, selection)
    if log(rand()) < weight
        # accept
        return new_trace
    else
        # reject
        return trace
    end
end

function rjmcmc(model, forward, forward_args_rest, backward, backward_args_rest,
                injective, injective_args, trace, correction)
    model_args = get_call_record(trace).args
    model_score = get_call_record(trace).score
    forward_args = (trace, forward_args_rest...,)
    forward_trace = simulate(forward, forward_args)
    forward_score = get_call_record(forward_trace).score
    input = pair(get_assignment(trace), get_assignment(forward_trace), :model, :proposal)
    (output, logabsdet) = apply(injective, injective_args, input)
    (model_constraints, backward_constraints) = unpair(output, :model, :proposal)
    new_trace = assess(model, model_args, model_constraints)
    new_model_score = get_call_record(new_trace).score
    backward_args = (new_trace, backward_args_rest...,)
    backward_trace = assess(backward, backward_args, backward_constraints)
    backward_score = get_call_record(backward_trace).score
    alpha = new_model_score - model_score - forward_score + backward_score + logabsdet + correction(new_trace)
    if log(rand()) < alpha
        # accept
        return new_trace
    else
        # reject
        return trace
    end
end

function sample_momenta(mass, n::Int)
    Float64[random(normal, 0, mass) for _=1:n]
end

function assess_momenta(momenta, mass)
    logprob = 0.
    for val in momenta
        logprob += logpdf(normal, val, 0, mass)
    end
    logprob
end

function hmc(model::Generator{T,U}, selection::AddressSet, trace::U;
             mass=0.1, L=10, eps=0.1) where {T,U}
    prev_model_score = get_call_record(trace).score
    model_args = get_call_record(trace).args

    # run leapfrog dynamics
    new_trace = trace
    local prev_momenta_score::Float64
    local momenta::Vector{Float64}
    for step=1:L

        # half step on momenta
        (_, values_trie, gradient_trie) = backprop_trace(model, new_trace, selection, nothing)
        values = to_array(values_trie, Float64)
        gradient = to_array(gradient_trie, Float64)
        if step == 1
            momenta = sample_momenta(mass, length(values))
            prev_momenta_score = assess_momenta(momenta, mass)
        else
            momenta += (eps / 2) * gradient
        end

        # full step on positions
        values_trie = from_array(values_trie, values + eps * momenta)

        # half step on momenta
        (new_trace, _, _) = update(model, model_args, NoChange(), new_trace, values_trie)
        (_, _, gradient_trie) = backprop_trace(model, new_trace, selection, nothing)
        gradient = to_array(gradient_trie, Float64)
        momenta += (eps / 2) * gradient
    end

    # assess new model score (negative potential energy)
    new_model_score = get_call_record(new_trace).score

    # assess new momenta score (negative kinetic energy)
    new_momenta_score = assess_momenta(-momenta, mass)

    # accept or reject
    alpha = new_model_score - prev_model_score + new_momenta_score - prev_momenta_score
    if log(rand()) < alpha
        new_trace
    else
        trace
    end
end

function mala(model::Generator{T,U}, selection::AddressSet, trace::U, tau) where {T,U}
    model_args = get_call_record(trace).args
    std = sqrt(2 * tau)

    # forward proposal
    (_, values_trie, gradient_trie) = backprop_trace(model, trace, selection, nothing)
    values = to_array(values_trie, Float64)
    gradient = to_array(gradient_trie, Float64)
    forward_mu = values + tau * gradient
    forward_score = 0.
    proposed_values = Vector{Float64}(undef, length(values))
    for i=1:length(values)
        proposed_values[i] = random(normal, forward_mu[i], std)
        forward_score += logpdf(normal, proposed_values[i], forward_mu[i], std)
    end

    # evaluate model weight
    constraints = from_array(values_trie, proposed_values)
    (new_trace, weight, discard) = update(
        model, model_args, NoChange(), trace, constraints)

    # backward proposal
    (_, _, backward_gradient_trie) = backprop_trace(model, new_trace, selection, nothing)
    backward_gradient = to_array(backward_gradient_trie, Float64)
    @assert length(backward_gradient) == length(values)
    backward_score = 0.
    backward_mu  = proposed_values + tau * backward_gradient
    for i=1:length(values)
        backward_score += logpdf(normal, values[i], backward_mu[i], std)
    end

    # accept or reject
    alpha = weight - forward_score + backward_score
    if log(rand()) < alpha
        return new_trace
    else
        return trace
    end
end

export mh, rjmcmc, hmc, mala

#######################
# importance sampling #
#######################

function importance_sampling(model::Generator{T,U}, model_args::Tuple,
                             observations::Assignment,
                             num_samples::Int) where {T,U}
    traces = Vector{U}(undef, num_samples)
    log_weights = Vector{Float64}(undef, num_samples)
    for i=1:num_samples
        (traces[i], log_weights[i]) = generate(model, model_args, observations)
    end
    log_total_weight = logsumexp(log_weights)
    log_ml_estimate = log_total_weight - log(num_samples)
    log_normalized_weights = log_weights .- log_total_weight
    return (traces, log_normalized_weights, log_ml_estimate)
end

function importance_sampling(model::Generator{T,U}, model_args::Tuple,
                             observations::Assignment,
                             proposal::Generator, proposal_args::Tuple,
                             num_samples::Int) where {T,U}
    traces = Vector{U}(undef, num_samples)
    log_weights = Vector{Float64}(undef, num_samples)
    for i=1:num_samples
        proposal_trace = simulate(proposal, proposal_args)
        proposal_score = get_call_record(proposal_trace).score
        constraints = merge(observations, get_assignment(proposal_trace))
        traces[i] = assess(model, model_args, constraints)
        model_score = get_call_record(traces[i]).score
        log_weights[i] = model_score - proposal_score
    end
    log_total_weight = logsumexp(log_weights)
    log_ml_estimate = log_total_weight - log(num_samples)
    log_normalized_weights = log_weights - log_total_weight
    return (traces, log_normalized_weights, log_ml_estimate)
end

function importance_resampling(model::Generator{T,U}, model_args::Tuple,
                               observations::Assignment,
                               num_samples::Int; verbose=false)  where {T,U,V,W}
    (model_trace::U, log_weight) = generate(model, model_args, observations)
    log_total_weight = log_weight
    for i=2:num_samples
        verbose && println("sample: $i of $num_samples")
        (cand_model_trace, log_weight) = generate(model, model_args, observations)
        log_total_weight = logsumexp(log_total_weight, log_weight)
        if bernoulli(exp(log_weight - log_total_weight))
            model_trace = cand_model_trace
        end
    end
    log_ml_estimate = log_total_weight - log(num_samples)
    return (model_trace::U, log_ml_estimate::Float64)
end

function importance_resampling(model::Generator{T,U}, model_args::Tuple,
                               observations::Assignment,
                               proposal::Generator{V,W}, proposal_args::Tuple,
                               num_samples::Int; verbose=false)  where {T,U,V,W}
    proposal_trace::W = simulate(proposal, proposal_args)
    proposal_score = get_call_record(proposal_trace).score
    constraints::Assignment = merge(observations, get_assignment(proposal_trace))
    (model_trace::U, model_score) = generate(model, model_args, constraints)
    log_total_weight = model_score - proposal_score
    for i=2:num_samples
        verbose && println("sample: $i of $num_samples")
        proposal_trace = simulate(proposal, proposal_args)
        proposal_score = get_call_record(proposal_trace).score
        constraints = merge(observations, get_assignment(proposal_trace))
        (cand_model_trace, model_score) = generate(model, model_args, constraints)
        log_weight = model_score - proposal_score
        log_total_weight = logsumexp(log_total_weight, log_weight)
        if bernoulli(exp(log_weight - log_total_weight))
            model_trace = cand_model_trace
        end
    end
    log_ml_estimate = log_total_weight - log(num_samples)
    return (model_trace::U, log_ml_estimate::Float64)
end

export importance_sampling, importance_resampling

######################
# particle filtering #
######################

function effective_sample_size(log_weights::Vector{Float64})
    # assumes weights are normalized
    log_ess = -logsumexp(2. * log_weights)
    exp(log_ess)
end

"""
the first argument to model should be an integer, starting from 1, that indicates the step
get_observations is a function of the step that returns a choice trie
rejuvenation_move is a function of the step and the previous trace, that returns a new trace
"""
function particle_filter(model::Generator{T,U}, num_steps::Int,
                         num_particles::Int, ess_threshold::Real,
                         get_observations::Function,
                         rejuvenation_move::Function=(t,trace) -> trace) where {T,U}

    log_unnormalized_weights = Vector{Float64}(undef, num_particles)
    log_ml_estimate = 0.
    observations = get_observations(1)
    traces = Vector{T}(undef, num_particles)
    next_traces = Vector{T}(undef, num_particles)
    for i=1:num_particles
        (traces[i], log_unnormalized_weights[i]) = generate(model, (1,), observations)
    end
    for step=2:num_steps

        # rejuvenation moves
        for i=1:num_particles
            traces[i] = rejuvenation_move(step, traces[i])
        end

        # compute new weights
        log_total_weight = logsumexp(log_unnormalized_weights)
        log_normalized_weights = log_unnormalized_weights .- log_total_weight

        # resample
        if effective_sample_size(log_normalized_weights) < ess_threshold
            weights = exp.(log_normalized_weights)
            parents = rand(Distributions.Categorical(weights / sum(weights)), num_particles)
            log_ml_estimate += log_total_weight - log(num_particles)
            log_unnormalized_weights = zeros(num_particles)
        else
            parents = 1:num_particles
        end

        # extend by one time step
        observations = get_observations(step)
        args_change = nothing
        for i=1:num_particles
            parent = parents[i]
            (next_traces[i], weight) = extend(model, (step,), args_change, traces[i], observations)
            log_unnormalized_weights[i] += weight
        end
        tmp = traces
        traces = next_traces
        next_traces = tmp
    end

    # finalize estimate of log marginal likelihood
    log_total_weight = logsumexp(log_unnormalized_weights)
    log_normalized_weights = log_unnormalized_weights .- log_total_weight
    log_ml_estimate += log_total_weight - log(num_particles)
    return (traces, log_normalized_weights, log_ml_estimate)
end

function particle_filter(model::Generator{T,U}, model_args_rest::Tuple,
                         num_steps::Int, num_particles::Int, ess_threshold::Real,
                         get_init_observations_and_proposal_args::Function,
                         get_step_observations_and_proposal_args::Function,
                         init_proposal::Generator, step_proposal::Generator;
                         verbose::Bool=false) where {T,U}

    log_unnormalized_weights = Vector{Float64}(undef, num_particles)
    log_ml_estimate = 0.
    (observations, proposal_args) = get_init_observations_and_proposal_args()
    traces = Vector{U}(undef, num_particles)
    next_traces = Vector{U}(undef, num_particles)
    for i=1:num_particles
        proposal_trace = simulate(init_proposal, proposal_args)
        proposal_score = get_call_record(proposal_trace).score
        constraints = merge(observations, get_assignment(proposal_trace))
        (traces[i], model_weight) = generate(model, (1, model_args_rest...), constraints) # make this assess when debugging?
        log_unnormalized_weights[i] = model_weight - proposal_score
    end

    local parents::Vector{Int}
    for step=2:num_steps

        # compute new weights
        log_total_weight = logsumexp(log_unnormalized_weights)
        log_normalized_weights = log_unnormalized_weights .- log_total_weight
        ess = effective_sample_size(log_normalized_weights)
        verbose && println("step: $step, ess: $ess")

        # resample
        if ess < ess_threshold
            verbose && println("resampling..")
            weights = exp.(log_normalized_weights)
            parents = rand(Distributions.Categorical(weights / sum(weights)), num_particles)
            log_ml_estimate += log_total_weight - log(num_particles)
            log_unnormalized_weights = zeros(num_particles)
        else
            parents = collect(1:num_particles)
        end

        # extend by one time step
        args_change = nothing # TODO establish a convention for args change (only the step argument changes)
        for i=1:num_particles
            parent = parents[i]
            parent_trace = traces[parent]
            (observations, proposal_args) = get_step_observations_and_proposal_args(step, parent_trace)
            proposal_trace = simulate(step_proposal, proposal_args)
            proposal_score = get_call_record(proposal_trace).score
            constraints = merge(observations, get_assignment(proposal_trace))
            (next_traces[i], model_weight) = extend(
                model, (step, model_args_rest...), args_change, parent_trace, constraints)
            log_weight = model_weight - proposal_score
            log_unnormalized_weights[i] += log_weight
        end
        tmp = traces
        traces = next_traces
        next_traces = tmp
    end

    # finalize estimate of log marginal likelihood
    log_total_weight = logsumexp(log_unnormalized_weights)
    log_normalized_weights = log_unnormalized_weights .- log_total_weight
    ess = effective_sample_size(log_normalized_weights)
    log_ml_estimate += log_total_weight - log(num_particles)
    verbose && println("final ess: $ess, final log_ml_est: $log_ml_estimate")
    return (traces, log_normalized_weights, log_ml_estimate)
end

export particle_filter

##############################
# Maximum a posteriori (MAP) #
##############################

"""
Backtracking gradient ascent for MAP inference on selected real-valued choices
"""
function map_optimize(model::Generator, selection::AddressSet,
                      trace; max_step_size=0.1, tau=0.5, min_step_size=1e-16, verbose=false)
    model_args = get_call_record(trace).args
    (_, values, gradient) = backprop_trace(model, trace, selection, nothing)
    values_vec = to_array(values, Float64)
    gradient_vec = to_array(gradient, Float64)
    step_size = max_step_size
    score = get_call_record(trace).score
    while true
        new_values_vec = values_vec + gradient_vec * step_size
        values = from_array(values, new_values_vec)
        # TODO discard and weight are not actually needed, there should be a more specialized variant
        (new_trace, _, discard, _) = update(model, model_args, NoChange(), trace, values)
        new_score = get_call_record(new_trace).score
        change = new_score - score
        if verbose
            println("step_size: $step_size, prev score: $score, new score: $new_score, change: $change")
        end
        if change >= 0.
            # it got better, return it
            return new_trace
        elseif step_size < min_step_size
            # it got worse, but we ran out of attempts
            return trace
        end

        # try again with a smaller step size
        step_size = tau * step_size
    end
end

export map_optimize

################
# SGD training #
################

struct SGDTrainConf
    num_batch::Int
    batch_size::Int
    num_minibatch::Int
    minibatch_size::Int
    input_extractor::Function
    constraint_extractor::Function
    minibatch_callback::Function
    batch_callback::Function
end

function sgd_train_batch(teacher::Generator{T,U}, teacher_args::Tuple,
                         batch_student::Generator{V,W}, conf::SGDTrainConf,
                         verbose=false) where {T,U,V,W}

    for batch=1:conf.num_batch

        # generate training batch
        training_assignment = Vector{Any}(undef, conf.batch_size)
        for i=1:conf.batch_size
            training_assignment[i] = get_assignment(simulate(teacher, teacher_args))
            if verbose && (i % 100 == 0)
                println("batch $batch, generating training data $i of $(conf.batch_size))")
            end
        end

        # train on this batch
        for minibatch=1:conf.num_minibatch
            permuted = Random.randperm(conf.batch_size)
            minibatch_idx = permuted[1:conf.minibatch_size]
            assignments = training_assignment[minibatch_idx]
            input = conf.input_extractor(assignments)
            constraints = conf.constraint_extractor(assignments)
            student_trace = assess(batch_student, input, constraints)
            avg_score = get_call_record(student_trace).score / conf.minibatch_size
            backprop_params(batch_student, student_trace, nothing)
            conf.minibatch_callback(batch, minibatch, avg_score, verbose)
        end

        conf.batch_callback(batch, verbose)
    end
end

export sgd_train_batch, SGDTrainConf

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
"""
$(SIGNATURES)

Find `\\newcommand` elements and try to parse what follows to form a proper Latex command.
Return a list of such elements.

The format is:
    \\newcommand{NAMING}[NARG]{DEFINING}
where [NARG] is optional (see `LX_NARG_PAT`).
"""
function find_md_lxdefs(tokens::Vector{Token}, blocks::Vector{OCBlock})
    # container for the definitions
    lxdefs = Vector{LxDef}()
    # find braces `{` and `}`
    braces = filter(β -> β.name == :LXB, blocks)
    nbraces = length(braces)
    # keep track of active tokens
    active_tokens = ones(Bool, length(tokens))
    active_blocks = ones(Bool, length(blocks))

    # go over active tokens, stop over the ones that indicate a newcommand
    # deactivate the tokens that are within the scope of a newcommand
    for (i, τ) ∈ enumerate(tokens)
        # skip inactive tokens
        active_tokens[i] || continue
        # look for tokens that indicate a newcommand
        (τ.name == :LX_NEWCOMMAND) || continue

        # find first brace blocks after the newcommand (naming)
        fromτ = from(τ)
        k = findfirst(β -> (fromτ < from(β)), braces)
        # there must be two brace blocks after the newcommand (name, def)
        if isnothing(k) || !(1 <= k < nbraces)
            throw(LxDefError("Ill formed newcommand (needs two {...})"))
        end

        # try to find a number of arg between these two first {...} to see
        # if it may contain something which we'll try to interpret as [.d.]
        rge = (to(braces[k])+1):(from(braces[k+1])-1)
        lxnarg = 0
        # it found something between the naming brace and the def brace
        # check if it looks like [.d.] where d is a number and . are
        # optional spaces (specification of the number of arguments)
        if !isempty(rge)
            lxnarg = match(LX_NARG_PAT, subs(str(braces[k]), rge))
            if isnothing(lxnarg)
                throw(LxDefError("Ill formed newcommand (where I expected the "*
                                 "specification of the number of arguments)."))
            end
            matched = lxnarg.captures[2]
            lxnarg = isnothing(matched) ? 0 : parse(Int, matched)
        end

        # assign naming / def
        naming_braces = braces[k]
        defining_braces = braces[k+1]
        # try to find a valid command name in the first set of braces
        matched = match(LX_NAME_PAT, content(naming_braces))
        if isnothing(matched)
            throw(LxDefError("Invalid definition of a new command expected a command name " *
                             "of the form `\\command`."))
        end

        # keep track of the command name, definition and where it stops
        lxname = matched.captures[1]
        lxdef = content(defining_braces)
        todef = to(defining_braces)
        # post-process the def
        lxdef = ignore_starting_line_spaces(lxdef)
        # store the new latex command
        push!(lxdefs, LxDef(lxname, lxnarg, lxdef, fromτ, todef))

        # mark newcommand token as processed as well as the next token
        # which is necessarily the command name (brace token is inactive here)
        active_tokens[i] = false
        # mark any block (inc. braces!) starting in the scope as inactive
        for (i, isactive) ∈ enumerate(active_blocks)
            isactive || continue
            (fromτ ≤ from(blocks[i]) ≤ todef) && (active_blocks[i] = false)
        end
    end # of enumeration of tokens

    # filter out the stuff that's now marked as inactive by virtue of being
    # part of a newcommand definition (these things will be inspected later)
    tokens = tokens[active_tokens]
    blocks = blocks[active_blocks]
    # separate the braces from the rest of the blocks, they will be used
    # to define the lxcoms
    braces_mask = map(β -> β.name == :LXB, blocks)
    braces = blocks[braces_mask]
    blocks = blocks[@. ~braces_mask]

    return lxdefs, tokens, braces, blocks
end

"""
$(SIGNATURES)

Retrieve the reference pointing to a `LxDef` corresponding to a given `lxname`.
If no reference is found but `inmath=true`, we propagate and let KaTeX deal with it. If something
is found, the reference is returned and will be accessed further down.
"""
function retrieve_lxdefref(lxname::SubString, lxdefs::Vector{LxDef},
                           inmath::Bool=false, offset::Int=0)::Ref
    # find lxdefs with matching name
    ks = findall(δ -> (δ.name == lxname), lxdefs)
    # check that the def is before the usage
    fromlx = from(lxname) + offset
    filter!(k -> (fromlx > from(lxdefs[k])), ks)
    if isempty(ks)
        if !inmath
            throw(LxComError("Command '$lxname' was used before it was defined."))
        end
        # not found but inmath --> let KaTex deal with it
        return Ref(nothing)
    end
    return Ref(lxdefs, ks[end])
end

"""
$(SIGNATURES)

Find `\\command{arg1}{arg2}...` outside of `xblocks` and `lxdefs`.
"""
function find_md_lxcoms(tokens::Vector{Token}, lxdefs::Vector{LxDef},
                        braces::Vector{OCBlock}, offset::Int=0;
                        inmath::Bool=false)::Tuple{Vector{LxCom},Vector{Token}}
    # containers for the lxcoms
    lxcoms   = Vector{LxCom}()
    active_τ = ones(Bool, length(tokens))
    nbraces  = length(braces)

    # go over tokens, stop over the ones that indicate a command
    for (i, τ) ∈ enumerate(tokens)
        active_τ[i] || continue
        (τ.name == :LX_COMMAND) || continue
        # 1. look for the definition given the command name
        lxname   = τ.ss
        lxdefref = retrieve_lxdefref(lxname, lxdefs, inmath, offset)
        # will only be nothing in a 'inmath' --> no failure, just ignore token
        isnothing(lxdefref[]) && continue

        # 2. retrieve number of arguments
        lxnarg = getindex(lxdefref).narg
        # 2.a there are none
        if lxnarg == 0
            push!(lxcoms, LxCom(lxname, lxdefref))
            active_τ[i] = false

        # >> there is at least one argument --> find all of them
        else
            # spot where an opening brace is expected
            nxtidx = to(τ) + 1
            b1_idx = findfirst(β -> (from(β) == nxtidx), braces)
            # --> it needs to exist + there should be enough braces left
            if isnothing(b1_idx) || (b1_idx + lxnarg - 1 > nbraces)
                throw(LxComError("Command '$lxname' expects $lxnarg argument(s) and there " *
                                 "should be no space(s) between the command name and the first " *
                                 "brace: \\com{arg1}..."))
            end

            # --> examine candidate braces, there should be no spaces between
            #  braces to avoid ambiguities
            cand_braces = braces[b1_idx:b1_idx+lxnarg-1]
            for bidx ∈ 1:lxnarg-1
                if (to(cand_braces[bidx]) + 1 != from(cand_braces[bidx+1]))
                    throw(LxComError("Argument braces should not be separated by space(s): " *
                                     "\\com{arg1}{arg2}... Verify a '$lxname' command."))
                end
            end

            # all good, can push it
            from_c = from(τ)
            to_c   = to(cand_braces[end])
            str_c  = subs(str(τ), from_c, to_c)
            push!(lxcoms, LxCom(str_c, lxdefref, cand_braces))

            # deactivate tokens in the span of the command (will be
            # reparsed later)
            first_active = findfirst(τ -> (from(τ) > to_c), tokens[i+1:end])
            if isnothing(first_active)
                active_τ[i+1:end] .= false
            elseif first_active > 1
                active_τ[i+1:(i+first_active-1)] .= false
            end
        end
    end
    return lxcoms, tokens[active_τ]
end


text:
using .StandardLibrary, .BiochemicalLibrary

const DEFAULT_LIBRARY = merge(standard_library, biochemical_library)

addlibrary!(newlib) = merge!(DEFAULT_LIBRARY, newlib)

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
# Recreates the analysis from the original data files

using SMLMAssociationAnalysis_NCB
using FileIO
using StatsBase
using StatsPlots

### Load saved data

samplenames = ["A", "B", "C", "D"]

nexperiments = 2
nreplicates = 3
nsamples = 4
ncells = 10

outputdir = "output/optimize"
datapath = joinpath(outputdir, "results_optimize.jld2")

experimentresults = load(datapath)["experimentresults"]

### Following optimization requires running run_optimize.jl instead of run_original.jl

# robustness (local distance)
average_percentileranks = [[[[[mean(x) for x ∈ experimentresults[l][k][j][i].percentileranks_by_localdensity]
                                       for i ∈ eachindex(experimentresults[l][k][j])]
                                       for j ∈ eachindex(experimentresults[l][k])]
                                       for k ∈ eachindex(experimentresults[l])]
                                       for l ∈ eachindex(experimentresults)]

average_percentileranks1 = [[average_percentileranks[l][k][j][i][m] for i in 1:ncells
                                                                    for j in 1:nsamples
                                                                    for k in 1:nreplicates
                                                                    for l in 1:nexperiments]
                                                                    for m in 1:10]

mean_average_percentileranks = mean.(average_percentileranks1)
sem_average_percentileranks = sem.(average_percentileranks1)
ci95_average_percentileranks = 1.96 .* sem_average_percentileranks

xvals = [200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800, 2000]
import StatsPlots.mm
plot(xvals, average_percentileranks, line=(0.25))
plot!(xvals, mean_average_percentileranks - ci95_average_percentileranks, line=(stroke(10), :gray))
plot!(xvals, mean_average_percentileranks + ci95_average_percentileranks, line=(stroke(10), :gray))
plot!(xvals, mean_average_percentileranks, line=(stroke(10), :black), yaxis=(0:1), xaxis=(0:400:2000), tickfontsize=124, legend=:none, size=(1024,2048), left_margin=21mm, top_margin=6mm, right_margin=10mm, bottom_margin=4mm, xrotation=45, ywiden=false)

savefig(joinpath(outputdir, "choice of local area effect on percentile rank.png"))

# temporal cutoff (needs special data)

ch1_counts = [[[[[length(x) for x ∈ experimentresults[l][k][j][i].channels[1].molecules]
                            for i ∈ eachindex(experimentresults[l][k][j])]
                            for j ∈ eachindex(experimentresults[l][k])]
                            for k ∈ eachindex(experimentresults[l])]
                            for l ∈ eachindex(experimentresults)]

ch1_frac = [[[[[x ./ ch1_counts[l][k][j][i][1] for x ∈ ch1_counts[l][k][j][i]]
                                            for i ∈ eachindex(experimentresults[l][k][j])]
                                            for j ∈ eachindex(experimentresults[l][k])]
                                            for k ∈ eachindex(experimentresults[l])]
                                            for l ∈ eachindex(experimentresults)]

ch1_frac1 = [[ch1_frac[l][k][j][i][m] for i in 1:ncells
                                       for j in 1:nsamples
                                       for k in 1:nreplicates
                                       for l in 1:nexperiments]
                                       for m in 1:10]

mean_ch1_frac = mean.(ch1_frac1)
sem_ch1_frac = sem.(ch1_frac1)
ci95_ch1_frac = 1.96 .* sem_ch1_frac

xvals = [0, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]

plot(xvals, ch1_frac, line=(0.25))
plot!(xvals, mean_ch1_frac, line=(stroke(10), :black), yaxis=(0:1), xaxis=(0:1000:5000), tickfontsize=124, legend=:none, size=(1024,2048), left_margin=21mm, top_margin=6mm, right_margin=10mm, bottom_margin=4mm, ywiden=false, xrotation=45)
plot!(xvals, mean_ch1_frac + ci95_ch1_frac, line=(stroke(10), :gray))
plot!(xvals, mean_ch1_frac - ci95_ch1_frac, line=(stroke(10), :gray))

savefig(joinpath(outputdir, "choice of t_off selection on number of molecules 1.png"))

ch2_counts = [[[[[length(x) for x ∈ experimentresults[l][k][j][i].channels[2].molecules]
                            for i ∈ eachindex(experimentresults[l][k][j])]
                            for j ∈ eachindex(experimentresults[l][k])]
                            for k ∈ eachindex(experimentresults[l])]
                            for l ∈ eachindex(experimentresults)]

ch2_frac = [[[[[x ./ ch2_counts[l][k][j][i][1] for x ∈ ch2_counts[l][k][j][i]]
                                                for i ∈ eachindex(experimentresults[l][k][j])]
                                                for j ∈ eachindex(experimentresults[l][k])]
                                                for k ∈ eachindex(experimentresults[l])]
                                                for l ∈ eachindex(experimentresults)]

ch2_frac1 = [[ch2_frac[l][k][j][i][m] for i in 1:ncells
                                       for j in 1:nsamples
                                       for k in 1:nreplicates
                                       for l in 1:nexperiments]
                                       for m in 1:10]

mean_ch2_frac = mean.(ch2_frac1)
sem_ch2_frac = sem.(ch2_frac1)
ci95_ch2_frac = 1.96 .* sem_ch2_frac

xvals = [0, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]

plot(xvals, ch2_frac, line=(0.25))
plot!(xvals, mean_ch2_frac, line=(stroke(10), :black), yaxis=(0:1), xaxis=(0:1000:5000), tickfontsize=124, legend=:none, size=(1024,2048),  left_margin=21mm, top_margin=6mm, right_margin=10mm, bottom_margin=4mm, ywiden=false, xrotation=45)
plot!(xvals, mean_ch2_frac + ci95_ch2_frac, line=(stroke(10), :gray))
plot!(xvals, mean_ch2_frac - ci95_ch2_frac, line=(stroke(10), :gray))

savefig(joinpath(outputdir, "choice of t_off selection on number of molecules 2.png"))