Update graph

demo1.m

@@ -9,18 +9,21 @@
 D(2) = init_D_from_txt('solway1.txt');
 D(3) = init_D_from_txt('solway2.txt');
 D(4) = init_D_from_txt('schapiro.txt');
+D(5) = init_D_from_txt('lynn.txt');
 
-for i = 1:length(D)
+for i = 5:length(D)
+    tic 
     [samples, post] = sample(D(i), h);
     for j = 1:length(samples)
         H(i,j) = samples(j);
         P(i,j) = post(j);
     end
+    toc
 end
-%}
 
-%load four.mat;
-load four_repro.mat;
+save demo1.mat
+%}
+load demo1.mat;
 
 figure;
 

demo_generate_samples.m

@@ -0,0 +1,24 @@
+h.alpha = 1.5;
+
+D(1) = init_D_from_txt('hourglass1.txt');
+
+nsubjects = 40;
+
+count_4_1 = 0;
+count_4_7 = 0;
+both = 0;
+for i = 1:nsubjects
+    [~, samples, ~] = sample_graph_update(D(1), h);
+    H = samples(1);
+    if (H.c(4) == H.c(1))
+        count_4_1 = count_4_1 + 1;
+    end
+    if (H.c(4) == H.c(7))
+        count_4_7 = count_4_7 + 1;
+    end
+    if (H.c(4) == H.c(7)) && (H.c(4) == H.c(1))
+       both = both + 1; 
+    end
+end
+
+save count_expirement_1.mat

demo_update_graph.m

@@ -0,0 +1,38 @@
+clear all;
+% {
+h.alpha = 1.5;
+
+D(1) = init_D_from_txt('hourglass.txt');
+
+for i = 1:1
+    tic 
+    [D(i), samples, post] = sample_graph_update(D(i), h);
+    for j = 1:length(samples)
+        H(i,j) = samples(j);
+        P(i,j) = post(j);
+    end
+    toc
+end
+
+save demo_update_graph.mat
+% }
+
+load demo_update_graph.mat;
+
+figure;
+
+k = 1;
+for i = 1:length(D)
+    post = P(i,:);
+    [~,I] = maxk(post, k);
+    for j = 1:k
+        subplot(length(D),k, (i-1)*k+j);
+        plot_H(H(i,I(j)), D(i));
+        if j == ceil(k/2)
+            %ylabel(D(i).name);
+            title(D(i).name);
+        end
+        set(gca, 'xtick', []);
+        set(gca, 'ytick', []);
+    end
+end

hourglass.txt

@@ -1,10 +1,19 @@
 Hourglass
-6 7
-1 3
-3 5
-1 5
-3 4
-4 2
-2 6
-4 6
+6 15
+1 3 1
+3 5 1
+1 5 1
+3 4 1
+4 2 1
+1 2 0
+1 4 0
+2 5 0
+3 6 0
+4 5 0
+5 6 0
+2 6 1
+1 6 0
+2 3 0
+4 6 1
+0
 0

hourglass1.txt

@@ -0,0 +1,37 @@
+Hourglass
+7 21
+1 2 0
+1 3 0
+1 4 0
+1 5 0
+1 6 0
+1 7 0
+2 3 0
+2 4 0
+2 5 0
+2 6 0
+2 7 0
+3 4 0
+3 5 0
+3 6 0
+3 7 0
+4 5 0
+4 6 0
+4 7 0
+5 6 0
+5 7 0
+6 7 0
+0
+12
+1 2 1
+1 3 1
+1 4 1
+2 3 1
+2 4 1
+3 4 1
+4 5 1
+4 6 1
+4 7 1
+5 6 1
+5 7 1
+6 7 1

hourglass2.txt

@@ -0,0 +1,37 @@
+Hourglass
+7 21
+1 2 0
+1 3 0
+1 4 0
+1 5 0
+1 6 0
+1 7 0
+2 3 0
+2 4 0
+2 5 0
+2 6 0
+2 7 0
+3 4 0
+3 5 0
+3 6 0
+3 7 0
+4 5 0
+4 6 0
+4 7 0
+5 6 0
+5 7 0
+6 7 0
+0
+12
+4 5 1
+4 6 1
+4 7 1
+5 6 1
+5 7 1
+6 7 1
+1 2 1
+1 3 1
+1 4 1
+2 3 1
+2 4 1
+3 4 1

init_D_from_txt.m

@@ -7,11 +7,14 @@
     N = A(1); M = A(2);
     D.G.N = N;
     D.G.E = zeros(N, N); % TODO sparse?
+    D.G.I = zeros(N, N);
     for k = 1:M
-        A = freadline(f, '%d %d');
-        i = A(1); j = A(2);
-        D.G.E(i,j) = 1;
-        D.G.E(j,i) = 1;
+        A = freadline(f, '%d %d %d');
+        i = A(1); j = A(2); exists = A(3);
+        D.G.E(i,j) = exists;
+        D.G.E(j,i) = exists;
+        D.G.I(i,j) = 1;
+        D.G.I(j,i) = 1;
     end
 
     D.tasks.s = [];
@@ -23,6 +26,16 @@
         D.tasks.s = [D.tasks.s s];
         D.tasks.g = [D.tasks.g g];
     end
+
+    % updates
+    A = freadline(f, '%d');
+    num_updates = A(1);
+    D.updates = [];
+    for k = 1:num_updates
+        A = freadline(f, '%d %d %d');
+        i = A(1); j = A(2); exists = A(3);
+        D.updates = [D.updates; i j exists];
+    end
 
     fclose(f);
 end

init_H.m

@@ -18,21 +18,6 @@
 
     H.p = betarnd(1,1); % TODO const 
     H.q = betarnd(1,1); % TODO const 
-    H.tp = betarnd(1,1); % TODO const 
-
-    H.N = length(cnt);
-    H.hp = betarnd(1,1); % TODO const 
-    H.E = zeros(H.N, H.N); % TODO sparse ?
-    for k = 1:H.N
-        for l = 1:k-1
-            if rand < H.hp
-                H.E(k,l) = 1;
-                H.E(l,k) = 1;
-            end
-        end
-    end
-
-    % TODO bridges
 end
 
 

loglik.m

@@ -1,10 +1,13 @@
 function logp = loglik(H, D, h)
-    % P(D|H) = P(G,tasks|H) = P(tasks|G,H) P(G|H) 
+    % P(D|H) = P(G|H)  
     %
 
     logp = 0;
     for i = 1:D.G.N
         for j = 1:i-1
+            if D.G.I(i, j) == 0
+                continue;
+            end
             if H.c(i) == H.c(j)
                 if D.G.E(i,j)
                     logp = logp + log(H.p);
@@ -22,15 +25,5 @@
             % TODO bridges
         end
     end
-
-    for i = 1:length(D.tasks.s)
-        s = D.tasks.s(i);
-        logp = logp + log(1 / D.G.N);
-
-        g = D.tasks.g(i);
-        P = ones(1, D.G.N);
-        P(H.c ~= H.c(s)) = H.tp;
-        logp = logp + log(P(g)) - log(sum(P));
-    end
 end
 

logprior.m

@@ -17,20 +17,8 @@
     end
     assert(isequal(cnt, H.cnt));
 
-    logp = logp + log(betapdf(H.p,1,1)) + log(betapdf(H.q,1,1)) + log(betapdf(H.tp,1,1)); % TODO const
+    logp = logp + log(betapdf(H.p,1,1)) + log(betapdf(H.q,1,1)); % TODO const
 
-    for k = 1:H.N
-        for l = 1:k-1
-            if H.E(k,l)
-                logp = logp + log(H.hp);
-            else
-                logp = logp + log(1 - H.hp);
-            end
-        end
-    end
-
-    logp = logp + log(betapdf(H.hp,1,1));
-
     % TODO bridges
 
 end

logsumexp.m

@@ -0,0 +1,19 @@
+function s = logsumexp(x,dim)
+
+    % Returns log(sum(exp(x),dim)) while avoiding numerical underflow.
+    % Default is dim = 1 (columns).
+
+    if nargin == 1
+        % Determine which dimension sum will use
+        dim = find(size(x)~=1,1);
+        if isempty(dim), dim = 1; end
+    end
+
+    % subtract the largest in each column
+    y = max(x,[],dim);
+    x = bsxfun(@minus,x,y);
+    s = y + log(sum(exp(x),dim));
+    i = find(~isfinite(y));
+    if ~isempty(i)
+        s(i) = y(i);
+    end

lynn.txt

@@ -0,0 +1,33 @@
+Lynn 2018
+15 30
+1 2
+2 3
+3 4
+4 5
+5 6
+6 7
+7 8
+8 9
+9 10
+10 11
+11 12
+12 13
+13 14
+14 15
+15 1
+1 3
+2 4
+3 5
+4 6
+5 7
+6 8
+7 9
+8 10
+9 11
+10 12
+11 13
+12 14
+13 15
+14 1
+15 2
+0

read_samples.m

@@ -0,0 +1,10 @@
+load count_expirement_1.mat
+
+% 4 is connected to 5, 6, 7
+
+n = nsubjects;
+c1 = count_4_1;
+c2 = count_4_7;
+
+p = 2 * binocdf(min(c1,c2), n, 0.5);
+y = binoinv([0.025 0.975], n, 0.5);