inst/coding_exploration.R

In tom-christie/transmit: Transmitting information using Poisson processes and Bayesian inference

get_angle <- function(a,b){acos( sum(a*b) / ( sqrt(sum(a * a)) * sqrt(sum(b * b)) ) )}


results <- data_frame()
for(i in 1:10000){
    signal_power <- 100
    n_neurons <- 2
    noise_power <- 1
    n_samples <- 5000
    n_codewords <- 2
    # each codeword is mapped to a different allocation
    #each column is a neuron, each row is a codeword
    allocation <- matrix(runif(n_neurons*n_codewords), ncol=n_neurons)
    allocation <- allocation/rowSums(allocation) * signal_power #normalize
    allocation <- allocation + noise_power #add noise
    #get a single sample from each neuron
    allocation <- matrix(c(0, 10, 10, 0),ncol=n_neurons,byrow=T)
    samples <- matrix(sapply(X = allocation, FUN = function(lambda) rpois(n_samples, lambda = lambda)),ncol=n_neurons,byrow=FALSE)
    message_sent <- rep(1:n_codewords,each=n_samples)
    #we want to calculate the likelihood of each row (pair of examples)
    #each row was generated by a codeword
    #and we want to calculate the likelihood that each row was generated by each codeword
    #that'll end up being n_codewords*n_codewords
    #want to end up with a matrix, samples x probability generated by that codeword
    #likelihood of each observation given the first rate
    #likelihood of each observation given the second rate
    likelihoods <- sapply(1:n_codewords,function(x) rowSums(dpois(samples, lambda = allocation[x,])))
    likelihoods <- likelihoods/rowSums(likelihoods)
    colnames(likelihoods) <- paste0('received ',1:n_codewords)
    rownames(likelihoods) <- paste0('sent ',message_sent)
    message_decoded <- unname(apply(likelihoods, MARGIN=1, FUN=which.max))
    correctly_decoded <- data_frame(
        message_sent = message_sent,
        message_decoded = message_decoded
    ) %>% 
        mutate(correct = (message_sent == message_decoded)) %>%
        count(correct) %>%
        filter(correct) %>%
        pull(n)
    accuracy = correctly_decoded/(n_samples*n_codewords)
    # if(accuracy > 0.9){
    #     print('good decoding')
    #     print(allocation)
    # }
    # if(get_angle(allocation[1,],allocation[2,]) > 1.25){
    #     print('high angle')
    #     print(allocation)
    # }
    results <- results %>% bind_rows(data_frame(
            accuracy = accuracy,
            angle = get_angle(allocation[1,],allocation[2,]),
            distance = dist(allocation),
            c1n1 = allocation[1,1],
            c1n2 = allocation[1,2],
            c2n1 = allocation[2,1],
            c2n2 = allocation[2,2]
        ))
}
results %>%
    ggplot() + 
    geom_point(aes(c1n1-c1n2,c2n1-c2n2,color=accuracy)) 


# sapply(c(1,2), 
#        function(index){
#            rowSums(cbind(likelihood1[,index], noise_log_likelihood[,-index]))
#        })