refreshRoutingMatrix.m

function [rt, rtfun, rtnodes, sn] = refreshRoutingMatrix(self, rates)
% [RT, RTFUN, CSMASK, RTNODES, SN] = REFRESHROUTINGMATRIX(RATES)
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
sn = self.sn;
if nargin == 1
    if isempty(sn)
        line_error(mfilename,'refreshRoutingMatrix cannot retrieve station rates, pass them as an input parameters.');
    else
        rates = sn.rates;
    end
end
M = sn.nnodes;
K = sn.nclasses;
arvRates = zeros(1,K);
stateful = find(sn.isstateful)';
 
indSource = find(sn.nodetype == NodeType.Source);
indOpenClasses = find(sn.njobs == Inf);
for r = indOpenClasses
    arvRates(r) = rates(sn.nodeToStation(indSource),r);
end
 
[rt, rtnodes, linksmat, chains] = self.getRoutingMatrix(arvRates);
sn = self.sn;
sn.chains = chains;
 
if self.enableChecks
    for r=1:K
        if all(sn.routing(:,r) == -1)
            line_error(mfilename,sprintf('Routing strategy in class %d is unspecified at all nodes.',r));
        end
    end
end
 
isStateDep = any(sn.isstatedep(:,3));
 
rnodefuncell = cell(M*K,M*K);
 
if isStateDep
    for ind=1:M % from
        for jnd=1:M % to
            for r=1:K
                for s=1:K
                    if sn.isstatedep(ind,3)
                        switch sn.routing(ind,r)
                            case RoutingStrategy.RROBIN
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_rr(ind, jnd, r, s, linksmat, state_before, state_after);
                            case RoutingStrategy.WRROBIN
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_wrr(ind, jnd, r, s, linksmat, state_before, state_after);
                            case RoutingStrategy.JSQ
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_jsq(ind, jnd, r, s, linksmat, state_before, state_after);
                            case RoutingStrategy.KCHOICES
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_kchoices(ind, jnd, r, s, linksmat, state_before, state_after);
                            case RoutingStrategy.RL
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_rl(ind, jnd, r, s, linksmat, state_before, state_after);
                            otherwise
                                rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(~,~) rtnodes((ind-1)*K+r, (jnd-1)*K+s);
                        end
                    else
                        rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(~,~) rtnodes((ind-1)*K+r, (jnd-1)*K+s);
                    end
                end
            end
        end
    end
end
 
statefulNodesClasses = [];
for ind=getIndexStatefulNodes(self)
    statefulNodesClasses(end+1:end+K)= ((ind-1)*K+1):(ind*K);
end
 
% we now generate the node routing matrix for the given state and then
% lump the states for non-stateful nodes so that run gives the routing
% table for stateful nodes only
statefulNodesClasses = [];
for ind=stateful
    statefulNodesClasses(end+1:end+K)= ((ind-1)*K+1):(ind*K);
end
 
if isStateDep
    rtfunraw = @(state_before, state_after) dtmc_stochcomp(cell2mat(cellfun(@(f) f(state_before, state_after), rnodefuncell,'UniformOutput',false)), statefulNodesClasses);
    rtfun = rtfunraw;
    %rtfun = memoize(rtfunraw); % memoize to reduce the number of stoch comp calls
    %rtfun.CacheSize = 6000^2;
else
    rtfun = @(state_before, state_after) dtmc_stochcomp(rtnodes, statefulNodesClasses);
end
 
nchains = size(chains,1);
inchain = cell(1,nchains);
for c=1:nchains
    inchain{c} = find(chains(c,:));
end
 
sn.rt = rt;
sn.rtnodes = rtnodes;
%sn.rtorig = RoutingMatrix.rtnodes2rtorig(sn); %% causes issues to
%JLINE.from_line_links in convering example_closedModel_6
sn.rtfun = rtfun;
sn.chains = chains;
sn.nchains = nchains;
sn.inchain = inchain;
for c=1:sn.nchains
    if range(sn.refstat(inchain{c}))>0
        line_error(mfilename,sprintf('Classes within chain %d (classes: %s) have different reference stations.',c,mat2str(find(sn.chains(c,:)))));
    end
end
self.sn = sn;
 
    function p = sub_rr(ind, jnd, r, s, linksmat, state_before, state_after)
        % P = SUB_RR(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER)
 
        isf = self.sn.nodeToStateful(ind);
        if isempty(state_before{isf})
            p = min(linksmat(ind,jnd),1);
        else
            if r==s
                p = double(state_after{isf}(sub_routeslot(ind, r, state_after{isf}))==jnd);
            else
                p = 0;
            end
        end
    end
 
    function slot = sub_routeslot(ind, r, state_i)
        % SLOT = SUB_ROUTESLOT(IND, R, STATE_I)
        % Column of STATE_I holding the round-robin pointer of class R at node
        % IND.
        %
        % The per-class routing pointers are appended in class order, and only
        % for the classes that actually route round-robin, followed by the
        % node-level block (nvars column 2R+1) if the node has one. The pointer
        % is therefore located by counting back from the end of the state: over
        % the node-level block, then over the pointers of the classes after r.
        %
        % Reading it as state_i(end-R+r) instead, as this did, assumes that the
        % pointers are exactly the last R columns. That silently reads the wrong
        % column whenever the node carries a trailing block of its own (the
        % polling controller of a POLLING station), and also whenever only some
        % of the classes route round-robin, in which case a pointer is confused
        % with a service phase.
        R = self.sn.nclasses;
        nodeblock = self.sn.nvars(ind, 2*R+1);
        after = sum(self.sn.nvars(ind, (R+r+1):(2*R)));
        slot = numel(state_i) - nodeblock - after;
    end
 
    function p = sub_wrr(ind, jnd, r, s, linksmat, state_before, state_after)
        % P = SUB_WRR(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER)
        % WRR slot holds a POSITION in weighted_outlinks; map it to the
        % destination node index before comparing.
 
        isf = self.sn.nodeToStateful(ind);
        if isempty(state_before{isf})
            p = min(linksmat(ind,jnd),1);
        else
            if r==s
                pos = state_after{isf}(sub_routeslot(ind, r, state_after{isf}));
                np = self.sn.nodeparam{ind}{r};
                if isfield(np,'weighted_outlinks') && ~isempty(np.weighted_outlinks) ...
                        && pos >= 1 && pos <= length(np.weighted_outlinks)
                    dest = np.weighted_outlinks(pos);
                elseif pos >= 1 && pos <= length(np.outlinks)
                    dest = np.outlinks(pos);
                else
                    p = 0; return;
                end
                p = double(dest == jnd);
            else
                p = 0;
            end
        end
    end
 
    function p = sub_jsq(ind, jnd, r, s, linksmat, state_before, state_after) %#ok<INUSD>
        % P = SUB_JSQ(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER) %#OK<INUSD>
 
        isf = self.sn.nodeToStateful(ind);
        if isempty(state_before{isf})
            p = min(linksmat(ind,jnd),1);
        else
            if r==s
                n = Inf*ones(1,self.sn.nnodes);
                for knd=1:self.sn.nnodes
                    if linksmat(ind,knd)
                        ksf = self.sn.nodeToStateful(knd);
                        n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                    end
                end
                if n(jnd) == min(n)
                    p = 1 / sum(n == min(n));
                else
                    p = 0;
                end
            else
                p = 0;
            end
        end
    end
 
 
    function p = sub_rl(ind, jnd, r, s, linksmat, state_before, state_after) %#ok<INUSD>
        % P = SUB_RL(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER) %#OK<INUSD>
 
        isf = self.sn.nodeToStateful(ind);
        if isempty(state_before{isf})
            p = min(linksmat(ind,jnd),1);
        else
            if r==s
                % ----- new added contents ----- %
                if self.nodes{ind}.output.outputStrategy{1,r}{5}==0         % state_size=0, use tabular value fn
                    % Multi-class: each class supplies its own value function.
                    % State.toMarginal aggregates queue lengths across classes,
                    % so the per-class call works as long as outputStrategy{r}
                    % is configured per class.
                    value_function = self.nodes{ind}.output.outputStrategy{1,r}{3};
                    nodes_need_action = self.nodes{ind}.output.outputStrategy{1,r}{4};
                    
                    if ~isempty(find(nodes_need_action==ind, 1))
                        indQueue = find(self.sn.nodetype == NodeType.Queue);
                        v = Inf*ones(1, self.sn.nnodes);                    % value fn
                        n = Inf*ones(1, self.sn.nnodes);                    % queue length
                        x = zeros(1, length(indQueue));                     % current state
                        for knd_idx=1:length(indQueue)
                            knd = indQueue(knd_idx);
                            ksf = self.sn.nodeToStateful(knd);              %% does the state removes the job from the departure node already?
                            x(knd_idx) = State.toMarginal(self.sn, knd, state_before{ksf});
                        end     
 
                        for knd = 1:self.sn.nnodes
                            if linksmat(ind, knd)
                                tmp = x+1;
                                tmp(indQueue == knd) = tmp(indQueue == knd) + 1;
                                if max(tmp) <= size(value_function, 1)
                                    ttmp = num2cell(tmp);
                                    v(knd) = value_function(ttmp{:});
                                end
                                ksf = self.sn.nodeToStateful(knd);
                                n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                            end
                        end
 
                        if min(v) < Inf && max(x+1) < size(value_function, 1)  % in action space
                            if v(jnd) == min(v)
                                p = 1 / sum(v == min(v));
                            else
                                p = 0;
                            end
                        else                                                % not in action space, use JSQ
                            if n(jnd) == min(n)
                                p = 1 / sum(n == min(n));
                            else
                                p = 0;
                            end
                        end
 
                    else                                                    % not in nodes_need_action: this node doesn't use RL results, use JSQ
                        n = Inf*ones(1,self.sn.nnodes);
                        for knd=1:self.sn.nnodes
                            if linksmat(ind,knd)
                                ksf = self.sn.nodeToStateful(knd);
                                n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                            end
                        end
                        if n(jnd) == min(n)
                            p = 1 / sum(n == min(n));
                        else
                            p = 0;
                        end
                    end
                
                elseif self.nodes{ind}.output.outputStrategy{1,r}{5}>0      % state_size>0, use fn approx for value fn
                    % Multi-class: class-specific coefficient row vector.
                    coeff = self.nodes{ind}.output.outputStrategy{1,r}{3};
                    nodes_need_action = self.nodes{ind}.output.outputStrategy{1,r}{4};
                    stateSize = self.nodes{ind}.output.outputStrategy{1,r}{5};
 
                    if ~isempty(find(nodes_need_action==ind, 1))
                        indQueue = find(self.sn.nodetype == NodeType.Queue);
                        v = Inf*ones(1, self.sn.nnodes);                    % value fn
                        n = Inf*ones(1, self.sn.nnodes);                    % queue length
                        x = zeros(1, length(indQueue));                     % current state
                        for knd_idx=1:length(indQueue)
                            knd = indQueue(knd_idx);
                            ksf = self.sn.nodeToStateful(knd);              %% does the state removes the job from the departure node already?
                            x(knd_idx) = State.toMarginal(self.sn, knd, state_before{ksf});
                        end     
                        for knd = 1:self.sn.nnodes
                            if linksmat(ind, knd)
                                tmp = x;
                                tmp(indQueue == knd) = tmp(indQueue == knd) + 1;
                                
                                tmp_vec = [1 tmp];
                                for i = 1:length(tmp)
                                    for j = i:length(tmp)
                                        tmp_vec(end+1) = tmp(i)* tmp(j);
                                    end
                                end
                                v(knd) = tmp_vec * coeff.';
 
                                ksf = self.sn.nodeToStateful(knd);
                                n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                            end
                        end
                        if min(v) < Inf  && max(x+1) < stateSize            % in action space
                            if v(jnd) == min(v)
                                p = 1 / sum(v == min(v));
                            else
                                p = 0;
                            end
                        else                                                % not in action space, use JSQ
                            if n(jnd) == min(n)
                                p = 1 / sum(n == min(n));
                            else
                                p = 0;
                            end
                        end
                    else                                                    % not in nodes_need_action: this node doesn't use RL results, use JSQ
                        n = Inf*ones(1,self.sn.nnodes);
                        for knd=1:self.sn.nnodes
                            if linksmat(ind,knd)
                                ksf = self.sn.nodeToStateful(knd);
                                n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                            end
                        end
                        if n(jnd) == min(n)
                            p = 1 / sum(n == min(n));
                        else
                            p = 0;
                        end     
                    end
 
                else                                                        % no value fn, use JSQ 
                % ----- end of new added contents ----- %
 
                    n = Inf*ones(1,self.sn.nnodes);
                    for knd=1:self.sn.nnodes
                        if linksmat(ind,knd)
                            ksf = self.sn.nodeToStateful(knd);
                            n(knd) = State.toMarginal(self.sn, knd, state_before{ksf});
                        end
                    end
                    if n(jnd) == min(n)
                        p = 1 / sum(n == min(n));
                    else
                        p = 0;
                    end
                end
            else
                p = 0;
            end
        end
    end
 
    function p = sub_kchoices(ind, jnd, r, s, linksmat, state_before, state_after) %#ok<INUSD>
        % P = SUB_KCHOICES(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER)
        % Power-of-K choices: matches LDES semantics in
        % Solver_ssj.kt:selectKChoicesDestinationWithClass.
        %
        % Without memory: enumerate the m^k ordered tuples obtained by sampling
        % k destinations WITH replacement, break ties by first occurrence in
        % the tuple, and return the marginal probability that jnd is chosen.
        %
        % With memory: the previously selected destination is forced as the
        % LAST candidate (LDES order), and the remaining k-1 candidates are
        % sampled iid with replacement. Memory is read from state_before; if
        % it is -1 (no prior pick) or no longer eligible, the closure falls
        % back to the memory-less form.
        %
        % Memoization: the marginal-probability vector p_dest(jnd) depends only
        % on (n[1..m], k, memPos). Across (ind, r) we cache by the per-call
        % queue-length signature so each unique (n, k, mem) is enumerated once
        % and the per-jnd probabilities are reused.
        persistent kchoicesCache
        if isempty(kchoicesCache)
            kchoicesCache = containers.Map('KeyType','char','ValueType','any');
        end
 
        isf = self.sn.nodeToStateful(ind);
        if isempty(state_before{isf})
            p = min(linksmat(ind,jnd),1);
            return;
        end
        if r ~= s
            p = 0; return;
        end
        eligible = find(linksmat(ind,:));
        m = numel(eligible);
        if m == 0 || ~any(eligible == jnd)
            p = 0; return;
        end
        np = self.sn.nodeparam{ind}{r};
        if isfield(np,'k') && ~isempty(np.k)
            k = np.k;
        else
            k = 2;
        end
        k = max(1, min(k, m));
        n = zeros(1, m);
        for i = 1:m
            knd = eligible(i);
            ksf = self.sn.nodeToStateful(knd);
            n(i) = State.toMarginal(self.sn, knd, state_before{ksf});
        end
        jnd_pos = find(eligible == jnd, 1);
 
        % Read memorized destination from the Router's nvars slot. The slot
        % stores the previously chosen destination as a node index (-1 = none).
        memEnabled = isfield(np, 'withMemory') && np.withMemory;
        memNode = -1;
        if memEnabled
            R = self.sn.nclasses;
            sb = state_before{isf};
            if numel(sb) >= R
                memNode = sb(end - R + r);
            end
        end
        memPos = -1;
        if memEnabled && memNode > 0
            tmpPos = find(eligible == memNode, 1);
            if ~isempty(tmpPos)
                memPos = tmpPos;
            end
        end
 
        % Cache key: (per-class) k, m, n vector, memPos. Same vector reused
        % across all jnd values for this (ind, r) pair.
        cacheKey = sprintf('k%d|m%d|mp%d|n%s', k, m, memPos, mat2str(n));
        if isKey(kchoicesCache, cacheKey)
            pVec = kchoicesCache(cacheKey);
        else
            pVec = zeros(1, m);
            if memPos > 0
                n_tuples = m^(k-1);
                tuple = zeros(1, k);
                tuple(end) = memPos;
                for ti = 0:(n_tuples-1)
                    rem = ti;
                    for c = 1:k-1
                        tuple(c) = mod(rem, m) + 1;
                        rem = floor(rem / m);
                    end
                    sub_n = n(tuple);
                    minval = min(sub_n);
                    winner_pos_in_tuple = find(sub_n == minval, 1);
                    winner = tuple(winner_pos_in_tuple);
                    pVec(winner) = pVec(winner) + 1;
                end
                pVec = pVec / n_tuples;
            else
                n_tuples = m^k;
                for ti = 0:(n_tuples-1)
                    tuple = zeros(1, k);
                    rem = ti;
                    for c = 1:k
                        tuple(c) = mod(rem, m) + 1;
                        rem = floor(rem / m);
                    end
                    sub_n = n(tuple);
                    minval = min(sub_n);
                    winner_pos_in_tuple = find(sub_n == minval, 1);
                    winner = tuple(winner_pos_in_tuple);
                    pVec(winner) = pVec(winner) + 1;
                end
                pVec = pVec / n_tuples;
            end
            % Bound cache size to avoid unbounded growth in long simulations.
            if kchoicesCache.Count > 50000
                remove(kchoicesCache, kchoicesCache.keys);
            end
            kchoicesCache(cacheKey) = pVec;
        end
        p = pVec(jnd_pos);
    end
 
end
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
% function [rt, rtfun, rtnodes, sn] = refreshRoutingMatrix(self, rates)
% % [RT, RTFUN, CSMASK, RTNODES, SN] = REFRESHROUTINGMATRIX(RATES)
% %
% % Copyright (c) 2012-2026, Imperial College London
% % All rights reserved.
% 
% sn = self.sn;
% if nargin == 1
%     if isempty(sn)
%         line_error(mfilename,'refreshRoutingMatrix cannot retrieve station rates, pass them as an input parameters.');
%     else
%         rates = sn.rates;
%     end
% end
% M = sn.nnodes;
% K = sn.nclasses;
% arvRates = zeros(1,K);
% stateful = find(sn.isstateful)';
% 
% indSource = find(sn.nodetype == NodeType.Source);
% indOpenClasses = find(sn.njobs == Inf);
% for r = indOpenClasses
%     arvRates(r) = rates(sn.nodeToStation(indSource),r);
% end
% 
% [rt, rtnodes, linksmat, chains] = self.getRoutingMatrix(arvRates);
% sn = self.sn;
% sn.chains = chains;
% 
% if self.enableChecks
%     for r=1:K
%         if all(sn.routing(:,r) == -1)
%             line_error(mfilename,sprintf('Routing strategy in class %d is unspecified at all nodes.',r));
%         end
%     end
% end
% 
% isStateDep = any(sn.isstatedep(:,3));
% 
% rnodefuncell = cell(M*K,M*K);
% 
% if isStateDep
%     for ind=1:M % from
%         for jnd=1:M % to
%             for r=1:K
%                 for s=1:K
%                     if sn.isstatedep(ind,3)
%                         switch sn.routing(ind,r)
%                             case RoutingStrategy.RROBIN
%                                 rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_rr(ind, jnd, r, s, linksmat, state_before, state_after);
%                             case RoutingStrategy.WRROBIN
%                                 rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_wrr(ind, jnd, r, s, linksmat, state_before, state_after);
%                             case RoutingStrategy.JSQ
%                                 rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(state_before, state_after) sub_jsq(ind, jnd, r, s, linksmat, state_before, state_after);
%                             otherwise
%                                 rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(~,~) rtnodes((ind-1)*K+r, (jnd-1)*K+s);
%                         end
%                     else
%                         rnodefuncell{(ind-1)*K+r, (jnd-1)*K+s} = @(~,~) rtnodes((ind-1)*K+r, (jnd-1)*K+s);
%                     end
%                 end
%             end
%         end
%     end
% end
% 
% statefulNodesClasses = [];
% for ind=getIndexStatefulNodes(self)
%     statefulNodesClasses(end+1:end+K)= ((ind-1)*K+1):(ind*K);
% end
% 
% % we now generate the node routing matrix for the given state and then
% % lump the states for non-stateful nodes so that run gives the routing
% % table for stateful nodes only
% statefulNodesClasses = [];
% for ind=stateful
%     statefulNodesClasses(end+1:end+K)= ((ind-1)*K+1):(ind*K);
% end
% 
% if isStateDep
%     rtfunraw = @(state_before, state_after) dtmc_stochcomp(cell2mat(cellfun(@(f) f(state_before, state_after), rnodefuncell,'UniformOutput',false)), statefulNodesClasses);
%     rtfun = rtfunraw;
%     %rtfun = memoize(rtfunraw); % memoize to reduce the number of stoch comp calls
%     %rtfun.CacheSize = 6000^2;
% else
%     rtfun = @(state_before, state_after) dtmc_stochcomp(rtnodes, statefulNodesClasses);
% end
% 
% nchains = size(chains,1);
% inchain = cell(1,nchains);
% for c=1:nchains
%     inchain{c} = find(chains(c,:));
% end
% 
% sn.rt = rt;
% sn.rtnodes = rtnodes;
% sn.rtfun = rtfun;
% sn.chains = chains;
% sn.nchains = nchains;
% sn.inchain = inchain;
% for c=1:sn.nchains
%     if range(sn.refstat(inchain{c}))>0
%         line_error(mfilename,sprintf('Classes within chain %d (classes: %s) have different reference stations.',c,mat2str(find(sn.chains(c,:)))));
%     end
% end
% self.sn = sn;
% 
%     function p = sub_rr(ind, jnd, r, s, linksmat, state_before, state_after)
%         % P = SUB_RR(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER)
% 
%         R = sn.nclasses;
%         isf = sn.nodeToStateful(ind);
%         if isempty(state_before{isf})
%             p = min(linksmat(ind,jnd),1);
%         else
%             if r==s
%                 p = double(state_after{isf}(end-R+r)==jnd);
%             else
%                 p = 0;
%             end
%         end
%     end
% 
%     function p = sub_wrr(ind, jnd, r, s, linksmat, state_before, state_after)
%         % P = SUB_WRR(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER)
% 
%         R = sn.nclasses;
%         isf = sn.nodeToStateful(ind);
%         if isempty(state_before{isf})
%             p = min(linksmat(ind,jnd),1);
%         else
%             if r==s
%                 p = double(state_after{isf}(end-R+r)==jnd);
%             else
%                 p = 0;
%             end
%         end
%     end
% 
%     function p = sub_jsq(ind, jnd, r, s, linksmat, state_before, state_after) %#ok<INUSD>
%         % P = SUB_JSQ(IND, JND, R, S, LINKSMAT, STATE_BEFORE, STATE_AFTER) %#OK<INUSD>
% 
%         isf = sn.nodeToStateful(ind);
%         if isempty(state_before{isf})
%             p = min(linksmat(ind,jnd),1);
%         else
%             if r==s
%                 n = Inf*ones(1,sn.nnodes);
%                 for knd=1:sn.nnodes
%                     if linksmat(ind,knd)
%                         ksf = sn.nodeToStateful(knd);
%                         n(knd) = State.toMarginal(sn, knd, state_before{ksf});
%                     end
%                 end
%                 if n(jnd) == min(n)
%                     p = 1 / sum(n == min(n));
%                 else
%                     p = 0;
%                 end
%             else
%                 p = 0;
%             end
%         end
%     end
% 
% end