some-julia

sims.jl

using Images
using Match

function sample_S(n, iters=1000)
    S = zeros(Int32, n)
    for step = 1:iters
        i = rand(1:n)
        @match (i,S[i]) begin
            (1,0),if S[i+1]==0 end              => S[i]=1
            (1,0),if S[i+1]==1 end              => continue
            (n,0),if S[i-1]==0 end              => S[i]=1
            (n,0),if S[i-1]==1 end              => continue
            (_,0),if S[i-1]==0 && S[i+1]==0 end => S[i]=1
            (_,1)                               => S[i]=0
            _                                   => continue
        end
    end
    S
end

function run_S(n=200, iters=1000000)
    P = pmap(sample_S,ones(Int32,iters)*n) # map w/ same param at each iter
    mean(P)
end

p = run_S()
# Note: Julia arrays are indexed by one => p[1]<=>p_0 and p[101]<=>p_100
println(size(p), ": p_0 = ", p[1], ", p_100 = ", p[101]) 
# (200,): p_0 = 0.333499, p_100 = 0.359054

function estimate_r(n)
    S = unique([sample_S(n) for i = 1:10000]) # probably enough iters to find all unique states for small values of n
    r = zeros(Int32,4)
    for i=1:length(S)
        @match S[i] begin
            [0,a...,0] => r[1] += 1
            [1,a...,0] => r[2] += 1
            [0,a...,1] => r[3] += 1
            [1,a...,1] => r[4] += 1
        end
    end
    r
end
[estimate_r(n) for n=2:20]
#19-element Array{Array{Int32,1},1}:
# Int32[1,1,1,0]            
# Int32[2,1,1,1]            
# Int32[3,2,2,1]            
# Int32[5,3,3,2]            
# Int32[8,5,5,3]            
# Int32[13,8,8,5]           
# Int32[21,13,13,8]         
# Int32[34,21,21,13]        
# Int32[55,34,34,21]        
# Int32[89,55,55,34]        
# Int32[144,89,89,55]       
# Int32[233,144,144,89]     
# Int32[377,233,233,144]    
# Int32[609,377,376,233]    
# Int32[964,601,594,371]    
# Int32[1433,894,890,551]   
# Int32[1992,1223,1227,763] 
# Int32[2481,1589,1518,947] 
# Int32[2911,1776,1818,1085]

function initial_checkerboard(n_particles=1278, size=100)
    board = zeros(size,size)
    
    # Randomly distribute the particles amongst the board
    for p = 1:n_particles
        while true
            i,j = rand(1:size), rand(1:size)
            if board[i,j] != 1
                board[i,j] = 1
                break
            end
        end
    end
    
    board
end

function run_checkerboard(iters=1000)
    S = initial_checkerboard()
    (_,n) = size(S)
    for step = 1:iters

        # each position
        for i=shuffle([1:n])
            for j=shuffle([1:n])
                if S[i,j]==1
                    direction = ["up", "left", "right", "down"][rand(1:4)]
                    @match direction begin
                        "up" => di, dj = i-1,j
                        "left" => di, dj = i,j+1
                        "right" => di, dj = i,j-1
                        "down" => di, dj = i+1,j
                    end
                    if di != 0 && di != n+1 && dj != 0 && dj != n+1 && S[di,dj] == 0
                        S[di,dj] = 1
                        S[i,j] = 0
                    end
                end
            end
        end
    end
    S
end

convert(Image, run_checkerboard())

function elemwise_mean(matrices, n, m)
    reshape(mean(map((i)-> reshape(i,n*m,1), matrices)), n,m)
    #reshape(mean(reshape(matrices',n*m,n_matrices),2), n,m)
end

A = [1 2; 3 4]
B = [5 6; 7 8]
C = [2 2; 2 2]
reshape(mean(reshape([A,B,C]',4,3),2), 2, 2)

A = [1 2; 3 4]
B = [5 6; 7 8]
C = [2 2; 2 2]

elemwise_mean([A,B,C], 2, 2)

P = [run_checkerboard() for i=1:1000]

#elemwise_mean(P, 1000, 100, 100)
#reshape(P', 10000, 1000)
avg_P = reshape(mean(map((m)-> reshape(m,10000,1), P)), 100,100)
convert(Image, avg_P)


N = 3
p = [0.5,0.5,0.5]

function neighbors(S,i)
    N = length(S)
    S[(i+1)%N + 1] + S[(i-1)%N + 1]
end

function sample_S(N,p)
    S = zeros(Int32, N)
    n_iter = 1000
    for step = 1:n_iter
        i = rand(1:N)
        m = neighbors(S,i)
        #print(S,' ',i,' ',m,'\n')
        if p[m+1] > rand()
            S[i] = 1
        else
            S[i] = 0
        end
    end
    S
end

function run(N, p, max_iterations=1000000)
    [sample_S(N,p) for i=1:max_iterations]
end

results = run(N, p)

counts = (String => Int64)[]

for r = results
    k = join(r)
    if haskey(counts,k)
        counts[k] += 1 
    else
        counts[k] = 0
    end
end

counts
#Dict{String,Int64} with 8 entries:
#  "000" => 125260
#  "111" => 124647
#  "001" => 125256
#  "011" => 125079
#  "101" => 123997
#  "110" => 125397
#  "010" => 125076
#  "100" => 125280