function [sv pv] = sbm_sample_v(theta, sh1)

    n = size(sh1, 1);
    raw = zeros(n, theta.num_visible, theta.L);
    
    repeated_a = repmat(reshape(theta.a, [1 size(theta.a)]), [n 1 1]);
    
    % contributions from layer 1 hiddens
    for p=1:theta.num_hidden1_patches

        stride_h1  = theta.num_hidden1 / theta.num_hidden1_patches;
        indices_h1 = (p-1)*stride_h1+1:p*stride_h1;
        indices_v  = theta.grid.patches(p, :);

        for l=1:theta.L
            raw(:, indices_v, l) = raw(:, indices_v, l) + sh1(:, indices_h1) * theta.W1(:, :, l);
        end

    end
    
    if theta.share_visible_biases
        
        % contributions from visible biases
        for p=1:theta.num_hidden1_patches

            stride_h1  = theta.num_hidden1 / theta.num_hidden1_patches;
            indices_v  = theta.grid.patches(p, :);

            raw(:, indices_v) = raw(:, indices_v) + repeated_a;

        end
        
    else
        
        raw = raw + repeated_a;
            
    end
        
    exp_raw = exp(raw);

    pv = exp_raw ./ repmat(sum(exp_raw, 3), [1 1 theta.L]);
    sv = reshape(draw_categorical(reshape(pv, [], theta.L)')', size(exp_raw));

end