In [ ]:

    

using Distributions


    

In [ ]:

    

function random_corr(num_traits::Int64, ke = 10.0^-3)
    random_corr = zeros(Float64, num_traits, num_traits)
    b = zeros(Float64, num_traits, num_traits)
    for i = 1:num_traits, j = 1:num_traits
        if(i >= j)
            b[i, j] = 1
        end
    end
    random_corr[2:num_traits, 1] = rand(Uniform(-1, 1), num_traits-1)
    b[2:num_traits, 1] = random_corr[2:num_traits, 1]
    for i in 2:num_traits
        b[i, 2:i] = sqrt(1 - random_corr[i, 1]^2)
    end

    for (i in 3:num_traits)
        for (j in 2:(i-1))
            b1 = dot(vec(b[j, 1:(j-1)]), vec(b[i, 1:(j-1)]))
            b2 = b[j, j] * b[i, j]
            z = b1 + b2
            y = b1 - b2
             if (b2 < ke)
                 random_corr[i, j] = b1
             else
                 random_corr[i, j] = y + (z - y) * rand(Uniform(0, 1))
             end
             cosinv = (random_corr[i, j] - b1)/b2
             if (isfinite(cosinv))
                 if (cosinv > 1)
                     b[i, (j+1):num_traits] = 0
                 elseif (cosinv < -1)
                     b[i, j] = -b[i, j]
                     b[i, (j+1):num_traits] = 0
                 else 
                     b[i, j] = b[i, j]*cosinv
                     sinTheta = sqrt(1 - cosinv^2)
                     b[i, (j+1):num_traits] = b[i, (j+1):num_traits]*sinTheta
                 end
             end
        end
    end
    random_corr = random_corr + random_corr' + diagm(ones(Float64, num_traits))
    perm = shuffle(collect(1:num_traits))
    random_corr[perm,:][:,perm]
end

function pcasimilarity(cov_x::Array{Float64, 2}, cov_y::Array{Float64, 2})
    eg_x = eig(cov_x)
    eg_y = eig(cov_y)
    eg_x_values = eg_x[1]
    eg_y_values = eg_y[1]
    eg_x_vectors = eg_x[2]
    eg_y_vectors = eg_y[2]

    total_var = dot(vec(eg_x_values), vec(eg_y_values))

    sum((eg_x_values * eg_y_values') .* (eg_x_vectors' * eg_y_vectors)^2)/total_var
end

function matrix_distance(cov_x, cov_y)
    sum(abs(cov_x - cov_y))
end


function mutation_b!(b_matrix, bin_b)
    mutation_b = rand(bin_b) + 1
    d_uni_p = DiscreteUniform(1, size(target_matrix, 1))
    d_uni_m = DiscreteUniform(1, size(target_matrix, 2))
    for k = range(1, mutation_b)
        i = rand(d_uni_p)
        j = rand(d_uni_m)
        b_matrix[i,j] = b_matrix[i, j] == 1 ? 0 : 1
    end
end

function project!(b_matrix, b_projection)
    copy!(b_projection, b_matrix * b_matrix')
    copy!(b_projection, b_projection ./ sqrt(diag(b_projection) * diag(b_projection)'))
end


    

In [ ]:

    

#target_matrix = diagm(5*ones(Float64, 10))
target_matrix  = abs(random_corr(5))


    

In [ ]:

    

b_2_dims = 100
p = size(target_matrix, 1)

t_init = 0.001
t_end = 0.0000001
t_steps = 1000000
n_monte = 100000 * p * b_2_dims

b_matrix = rand(Bernoulli(), size(target_matrix, 1), b_2_dims)

b_projection = ones(Float64, p, p)

d_uni_p = DiscreteUniform(1, size(target_matrix, 1))
d_uni_m = DiscreteUniform(1, size(target_matrix, 2))
uni_0_1 = Uniform(0, 1)


    

In [ ]:

    

while(matrix_distance(target_matrix, b_projection) > 2)
    for(t in linspace(t_init, t_end, t_steps))
        project!(b_matrix, b_projection)
        dist_before = matrix_distance(target_matrix, b_projection)

        i = rand(d_uni_p)
        j = rand(d_uni_m)
        b_matrix[i, j] = b_matrix[i, j] == 1 ? 0 : 1

        project!(b_matrix, b_projection)
        dist_after = matrix_distance(target_matrix, b_projection)

        if(dist_after < dist_before)
            #println(dist_after)
        elseif(exp((dist_before - dist_after)/t) < rand(uni_0_1))
            #println(dist_after)
        else
            b_matrix[i, j] = b_matrix[i, j] == 1 ? 0 : 1
        end
    end
end


    

In [ ]:

    

matrix_distance(target_matrix, b_projection)


    

In [ ]:

    

b_projection