function [sh1 ph1] = sbm_sample_h1(theta, sv, sh2, bottom_up_pass, compensation_factor)

    n = size(sv, 1);
    raw = zeros(n, theta.num_hidden1);
                           
    % contributions from visibles
    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_h1) = raw(:, indices_h1) + sv(:, indices_v, l) * theta.W1(:, :, l)';
        end

    end
        
    if nargin > 3 && bottom_up_pass
        
        % contributions from layer 1 biases
        raw = compensation_factor * raw + repmat_fast(theta.b1', n);
        
    else
       
        % contributions from layer 2 hiddens and layer 1 biases
        raw = raw + sh2 * theta.W2 + repmat_fast(theta.b1', n);
        
    end
    
    ph1 = sigmoid(raw/theta.T);
    sh1 = draw_bernoulli(ph1);
    
end