Imagining a 2-Hour, 4-Week Psalter

• Overview
• Comparing Major Psalters
  • Monastic (c. 540)
  • Pope St. Pius X (1911)
  • Liturgy of the Hours (1971)
• Analyzing Psalm Lengths
• Choosing Divisions
• Grouping Psalms into Hours
  • Goals and Approach
  • Implementation and First Attempt
  • Revisions and Final Result
• Filling the Calendar
  • Comparing Psalm Placement in Existing Psalters
• Analyzing the Results
• Next Steps and Future Work
• Notes and References

Overview

I have always been attracted to the idea of praying all 150 Psalms in one week. I have experimented with praying two such Psalters. First, the Anglican Breviary, which is simply an English translation of the 1911 Divino Afflatu Divine Office, using Pope St. Pius X’s 1-week, but very efficient Psalter. Later, the 1962 Monastic Diurnal (with supplemental resources for Matins), which still uses St. Benedict’s arrangement from the sixth century. Both proved to be too long and arduous for a simple and busy layman such as myself.

Most recently, I’ve transitioned to using the 4-volume Liturgy of the Hours. Despite disliking it more than both of the two previous versions of the Divine Office that I’ve used, it’s much shorter, and much more manageable. But due to some combination of busy-ness and a lack of discipline, I’ve only managed to be consistent with Lauds and Vespers, which leaves me with a 4-week cycle which does not even include the entire Psalter.

This all led to the question: Is a 4-week, 150-Psalm Psalter possible if we restrict ourselves to only Lauds and Vespers? Ideally, this arrangement would allow us to reuse the existing antiphons from the LotH, meaning that we’d only be able to pray 3 Psalms per hour.

Simple math says yes:

• 3 Psalms per hour * 2 hours per day =
• 6 Psalms per day * 7 days per week =
• 42 Psalms per week * 4 weeks per cycle =
• 168 Psalms

If we eliminate all of the canticles, it allows us to pray every Psalm, and gives us 18 slots to fill with divisions for the longer Psalms. This is many fewer division slots than the current LotH, so we will have to be strategic in how we use them.

Our goal is to create a 4-week, 2-hour Psalter that:

• Uses all 150 Psalms with minimal repetition (likely no repetition at all)
• Selects Psalm divisions that produce the greatest uniformity in length
• Arranges the Psalms in a way that conforms as much as possible to traditional Psalters

Comparing Major Psalters

Our first step is to review existing Psalters. For simplicity, I’ve chosen 3: the Monastic Psalter, Pope St. Pius X’s Psalter, and the Liturgy of the Hours’ Psalter. Hovering any of the Psalms in the tables highlights all of the instances of that Psalm in all the tables. Two notes about this:

1. Hovering highlights all instances of the Psalm, including divisions.
2. It is possible to switch between Greek and Masoretic numbering. But the highlighting is based on the Greek numbering. So, for example, hovering over an instance of Psalm 9 in the Masoretic numbering scheme will also highlight instances of Psalm 10, since Psalm 9 and 10 in the Masoretic numbering scheme are both the same Psalm in the Greek numbering scheme (Psalm 9).

Monastic (c. 540)

Pope St. Pius X (1911)

Liturgy of the Hours (1971)

Analyzing Psalm Lengths

The next step is to analyze the lengths of all the Psalms so that we can select which Psalms we want to divide, and how to divide them. For this (and for the rest of this project), I have used the Abbey Psalms and Canticles. Not because I have any attachment to or preference for this translation, but rather because this project was invisioned as an alternative Psalter for the Liturgy of the Hours, and starting in a year from the time this project began, this translation was going to begin being used with the new edition of the Liturgy of the Hours (2027).

Below are all the Psalms, sorted by the number of Verses.

We can look at some charts to help us visualize. First, the length of all the Psalms in order.

psalmNumbers = string(1:150);
numVerses = ...
    [ 6, 12,  9,  9, 13, 11, 18, 10, 21, 18,  7,  9,  6,   7,  5, ...
     11, 15, 51, 15, 10, 14, 32,  6, 10, 22, 12, 14,  9,  11, 13, ...
     25, 11, 22, 23, 28, 13, 40, 23, 14, 18, 14, 12,  5,  27, 18, ...
     12, 10, 15, 21, 23, 21, 11,  7,  9, 24, 14, 12, 12,  18, 14, ...
      9, 13, 12, 11, 14, 20,  8, 36, 37,  6, 24, 20, 28,  23, 11, ...
     13, 21, 72, 13, 20, 17,  8, 19, 13, 14, 17,  7, 19,  53, 17, ...
     16, 16,  5, 23, 11, 13, 12,  9,  9,  5,  8, 29, 22,  35, 45, ...
     48, 43, 14, 31,  7, 10, 10,  9,  8, 18, 19,  2, 29, 176,  7, ...
      8,  9,  4,  8,  5,  6,  5,  6,  8,  8,  3, 18,  3,   3, 21, ...
     26,  9,  8, 24, 14, 10,  8, 12, 15, 21, 10, 20, 14,   9,  6];
bar(psalmNumbers,numVerses,1,FaceColor='flat',EdgeColor='black',LineWidth=1);
xlabel('Psalm Number');
ylabel('Number of Verses');
title('Verses Per Psalm');

Next, we can sort them by length. Here, we add a bar (purple) for the average Psalm length. The average is ~16.8 verses.

avgPsalmLength = round(mean(numVerses));
numVersesWithAvg = [numVerses avgPsalmLength];
psalmNumbersWithAvg = [psalmNumbers 'average'];
[sortedNumVerses, sortIdx] = sort(numVersesWithAvg,'descend');
sortedPsalmNumbers = psalmNumbersWithAvg(sortIdx);
idxOfAvg = find(sortedPsalmNumbers=='average');
b = bar(sortedPsalmNumbers,sortedNumVerses,1,FaceColor='flat',EdgeColor='black',LineWidth=1);
b.CData(idxOfAvg,:) = [0.5 0 0.5];
xlabel('Psalm Number');
ylabel('Number of Verses');
title('Verses Per Psalm -- Sorted');

This is too many bars to be useful, so let’s zoom in on the 20 longest Psalms. The average is also included here, although there are many more than 20 Psalms longer than the average length.

n = 20;
topPsalmNumbers = sortedPsalmNumbers(1:n);
topNumVerses = sortedNumVerses(1:n);
topPsalmNumbersWithAvg = [topPsalmNumbers 'average'];
topNumVersesWithAvg = [topNumVerses avgPsalmLength];
b = bar(topPsalmNumbersWithAvg,topNumVersesWithAvg,1,FaceColor='flat',EdgeColor='black',LineWidth=1);
b.CData(end,:) = [0.5 0 0.5];
xlabel('Psalm Number');
ylabel('Number of Verses');
ylim([0 200]);
title(['Verses Per Psalm -- Sorted, Top ' num2str(n)]);
text(1:n+1,topNumVersesWithAvg,num2str(topNumVersesWithAvg'),'vert','bottom','horiz','center');

Choosing Divisions

Eighteen extra Psalm slots is not very much, so we want to use them as effectively as possible. It turns out that finding a good way to use them is not that hard. Some quick back-of-napkin math shows that after our divisions, our longest Psalm or Psalm division will be <35.

The first order of business is Psalm 119, which has 176 verses. This is a special Psalm, since it is a 22-part acrostic poem, with each part consisting of 8 verses. We do not want to split the sections, so we’ll have to break up the Psalm in chunks of 8. The Monastic Psalter and the Liturgy of the Hours breaks it into 22 divisions of 8 verses, and the Pope St. Pius X Psalter breaks it into 11 divisions of 16 verses. This is too many divisions for us, since we have rather limited real estate to work with. I opted for 7 divisions: 1 division with 32 verses and 6 with 24. This has the nice result of allowing us to pray the Psalm throughout an entire week. So for Psalm 119, we use an additional 6 Psalm slots.

The next Psalm is 78, with 72 verses. Two divisions of 36 verses are a little long, so we’ll divide it into 3 parts, using up 2 more Psalm slots.

Having used 8 slots, we have 10 left. The next longest Psalm is 89, which has 53 verses. So from here, we can split the next 9 longest Psalms into 2: 89, 18, 106, 105, 107, 37, 69, 68, 104. The next two longest Psalms are 22 (32 verses) and 109 (31 verses). Because Psalm 22 is such a beautiful and important Psalm, I’ve elected to keep it intact, and use the last Psalm slot to split Psalm 109.

The following table shows all of our divisions, and the resulting lengths. Note that the divisions are not always even. An effort was made to review the contents of the Psalm and divide them in a way that respects the structure and flow of the Psalm. Exactly how we divide them does not matter too much, since, as we will discuss later, we will want to force all of these divided Psalms to be prayed back-to-back in the same hour (with the exception of Psalm 119). So lopsided divisions will not affect the length of whole hours later on.

Below, we can see the original Top 20 Psalms, but now what they look like divided. As we can see, the average Psalm length is now ~15.0 verses.

topPsalmNumbersWithDivisions = ...
    ["119a" "119b" "119c" "119d" "119e" "119f" "119g" "78a" "78b" "78c" ...
     "89a" "89b" "18a" "18b" "106a" "106b" "105a" "105b" "107a" "107b" ...
     "37a" "37b" "69a" "69b" "68a" "68b" "104a" "104b" "22" "109a" ...
     "109b" "102" "118" "35" "73" "44" "136" "31" "average"];
topPsalmNumVersesWithDivisions = ...
    [32 24 24 24 24 24 24 31 24 17 19 34 25 26 23 25 22 23 22 ...
    21 22 18 18 19 19 17 23 12 32 20 11 29 29 28 28 27 26 25 15];
[sortedTopPsalmsNumVersesWithDivisions, sortIdx] = ...
    sort(topPsalmNumVersesWithDivisions,'descend');
sortedTopPsalmNumbersWithDivisions = topPsalmNumbersWithDivisions(sortIdx);
idxOfAvg = find(sortedTopPsalmNumbersWithDivisions=='average');
b = bar(sortedTopPsalmNumbersWithDivisions,sortedTopPsalmsNumVersesWithDivisions, ...
    1,FaceColor='flat',EdgeColor='black',LineWidth=1);
b.CData(idxOfAvg,:) = [0.5 0 0.5];
xlabel('Psalm Number');
ylabel('Number of Verses');
title(['Verses Per Psalm -- Sorted, Top 20, with Divisions']);
text(1:length(sortedTopPsalmsNumVersesWithDivisions),...
    sortedTopPsalmsNumVersesWithDivisions,...
    num2str(sortedTopPsalmsNumVersesWithDivisions'), ...
    'vert','bottom','horiz','center');

Now we can look at what the list of all Psalms now looks like. This next table shows all the Psalms, divided as specified, and sorted by length:

And finally, we can visualize all of these Psalms and divisions in one big chart:

allPsalmNumbers = ...
    ["89b" "119a" "22" "78a" "102" "118" "35" "73" "44" "18b" "136" ...
    "18a" "106b" "31" "119b" "119c" "119d" "119e" "119f" "119g" "78b" ...
    "55" "71" "139" "106a" "105b" "104a" "34" "38" "50" "74" "94" "105a" ...
    "107a" "37a" "25" "33" "103" "107b" "9" "49" "51" "77" "135" "145" ...
    "109a" "66" "72" "80" "147" "89a" "69b" "68a" "83" "88" "116" "37b" ...
    "69a" "7" "10" "40" "45" "59" "115" "132" "78c" "68b" "81" "86" "90" ...
    "91" "92" "average" "17" "19" "48" "144" "21" "27" "39" "41" "56" ...
    "60" "65" "85" "108" "140" "148" "5" "30" "36" "62" "76" "79" "84" ...
    "96" "104b" "2" "26" "42" "46" "57" "58" "63" "97" "143" "109b" "6" ...
    "16" "29" "32" "52" "64" "75" "95" "8" "20" "24" "47" "111" "112" ...
    "141" "146" "3" "4" "12" "28" "54" "61" "98" "99" "113" "122" "137" ...
    "149" "67" "82" "101" "114" "121" "124" "129" "130" "138" "142" "11" ...
    "14" "53" "87" "110" "120" "1" "13" "23" "70" "126" "128" "150" ...
    "15" "43" "93" "100" "125" "127" "123" "131" "133" "134" "117"];
allPsalmNumVerses = ...
    [34 32 32 31 29 29 28 28 27 26 26 25 25 25 24 24 24 24 24 24 24 24 ...
     24 24 23 23 23 23 23 23 23 23 22 22 22 22 22 22 21 21 21 21 21 21 ...
     21 20 20 20 20 20 19 19 19 19 19 19 18 18 18 18 18 18 18 18 18 17 ...
     17 17 17 17 16 16 15 15 15 15 15 14 14 14 14 14 14 14 14 14 14 14 ...
     13 13 13 13 13 13 13 13 12 12 12 12 12 12 12 12 12 12 11 11 11 11 ...
     11 11 11 11 11 10 10 10 10 10 10 10 10 9 9 9 9 9 9 9 9 9 9 9 9 8 8 ...
     8 8 8 8 8 8 8 8 7 7 7 7 7 7 6 6 6 6 6 6 6 5 5 5 5 5 5 4 3 3 3 2];
idxOfAvg = find(allPsalmNumbers=='average');
b = bar(allPsalmNumbers,allPsalmNumVerses,1, ...
    FaceColor='flat',EdgeColor='black',LineWidth=1);
b.CData(idxOfAvg,:) = [0.5 0 0.5];
xlabel('Psalm Number');
ylabel('Number of Verses');
title('Verses Per Psalm -- Sorted, with Divisions');

Grouping Psalms into Hours

At first, I had envisioned manually grouping Psalms into triples (one triple per hour). But as I thought about it, I realized that it would be tedious to get good results by hand.

At its core, this is an optimization problem. We have constraints and goals, and since we also have a way to measure how well a solution meets its goals, we can write a cost function.

Goals and Approach

To start out, we can enumerate some of our most constraints, those which define the problem:

1. Three Psalms are assigned to each group (or “triple”, or “bucket”).
2. Each Psalm is used exactly once.

We can impose additional constraints in order to form a solution to our liking:

1. All Psalm divisions must be grouped together, except:
2. Each division of Psalm 119 must be in its own group.
3. Triples should be no fewer than average-15 verses and no more than average+15 verses.

All the divided Psalms except for 119 have 3 or fewer parts, so constraint (3) can be accomplished. And we divided 119 into 7 groups, with the intention that each division be prayed on a different day for 7 days in a row.

In addition to constraints, we have some goals. For these, we very likely can’t accomplish them perfectly, but we can program our cost function so that we get the best solution possible. For our goals, we say that we are looking for a grouping of Psalms such that:

1. The length of each triple (total number of verses) is as close to the average as possible (~45 verses).
2. Wherever possible, we prefer buckets with sequential Psalms.

Finally, we can specify if there are some triples that we pre-define as requirements. To begin, I’ve selected a few, largely based on the Psalm groupings of the Monastic Psalter (and one triple which is reiteration of Constraint #3):

• 120, 121, 122
• 123, 124, 125
• 126, 127, 128
• 110, 111, 112
• 78a, 78b, 78c

Based on all of these constraints and goals, we can write a program which produces an optimal solution. We can then review the solution, and decide what changes we want to make, such as giving different weights to our goals, or loosening or tightening our constraints.

Implementation and First Attempt

Since I have a Computer Science background and not a math background, I had no idea how to begin programming for an optimization problem. Fortunately, ChatGPT was very helpful in pointing me in the right direction. Below is the first (correct) iteration of the implemention, which produces a good result. Here, I use MATLAB since I have access to all of its toolboxes, including the Optimization Toolbox.

N = 168;
psalm_names = ...
    ["89b" "119a" "22" "78a" "102" "118" "35" "73" "44" "18b" "136" ...
    "18a" "106b" "31" "119b" "119c" "119d" "119e" "119f" "119g" "78b" ...
    "55" "71" "139" "106a" "105b" "104a" "34" "38" "50" "74" "94" "105a" ...
    "107a" "37a" "25" "33" "103" "107b" "9" "49" "51" "77" "135" "145" ...
    "109a" "66" "72" "80" "147" "89a" "69b" "68a" "83" "88" "116" "37b" ...
    "69a" "7" "10" "40" "45" "59" "115" "132" "78c" "68b" "81" "86" "90" ...
    "91" "92" "17" "19" "48" "144" "21" "27" "39" "41" "56" ...
    "60" "65" "85" "108" "140" "148" "5" "30" "36" "62" "76" "79" "84" ...
    "96" "104b" "2" "26" "42" "46" "57" "58" "63" "97" "143" "109b" "6" ...
    "16" "29" "32" "52" "64" "75" "95" "8" "20" "24" "47" "111" "112" ...
    "141" "146" "3" "4" "12" "28" "54" "61" "98" "99" "113" "122" "137" ...
    "149" "67" "82" "101" "114" "121" "124" "129" "130" "138" "142" "11" ...
    "14" "53" "87" "110" "120" "1" "13" "23" "70" "126" "128" "150" ...
    "15" "43" "93" "100" "125" "127" "123" "131" "133" "134" "117"];
psalm_lengths = ...
    [34 32 32 31 29 29 28 28 27 26 26 25 25 25 24 24 24 24 24 24 24 24 ...
    24 24 23 23 23 23 23 23 23 23 22 22 22 22 22 22 21 21 21 21 21 21 ...
    21 20 20 20 20 20 19 19 19 19 19 19 18 18 18 18 18 18 18 18 18 17 ...
    17 17 17 17 16 16 15 15 15 15 14 14 14 14 14 14 14 14 14 14 14 ...
    13 13 13 13 13 13 13 13 12 12 12 12 12 12 12 12 12 12 11 11 11 11 ...
    11 11 11 11 11 10 10 10 10 10 10 10 10 9 9 9 9 9 9 9 9 9 9 9 9 8 8 ...
    8 8 8 8 8 8 8 8 7 7 7 7 7 7 6 6 6 6 6 6 6 5 5 5 5 5 5 4 3 3 3 2];

% These parameters can be tweaked to produce different results.
w_val = 1.0; % weight of value balancing objective
w_seq = 1.0; % weight of sequentiality objective
target_length = mean(psalm_lengths) * 3; % target length of a triple of Psalms
max_length = target_length + 15; % maximum length of a triple of Psalms
min_length = target_length - 15; % minimum length of a triple of Psalms

% Force Psalms split into pairs to go together, codified as a key-value
% pair. Allows us to ensure that if we see <key> in a triple, that the
% triple also contains <val>.
pairs = dictionary(...
    ["89a" "89b" "18a" "18b" "106a" "106b" "105a" "105b" "107a" "107b" ...
     "37a" "37b" "69a" "69b" "68a" "68b" "104a" "104b" "109a" "109b"], ...
     ["89b" "89a" "18b" "18a" "106b" "106a" "105b" "105a" "107b" "107a" ...
     "37b" "37a" "69b" "69a" "68b" "68a" "104b" "104a" "109b" "109a"]);

% Ensure that these triples of Psalms are grouped together.
forced_triples = [
    "120" "121" "122"
    "123" "124" "125"
    "126" "127" "128"
    "110" "111" "112"
    "78a" "78b" "78c"
];


%% Preprocessing
% The first step is to perform some simple transformations on the data so
% that it can more easily be processed.

% Exclude all "forced" Psalms from lists
all_indices = 1:N;
excluded_psalms = unique(forced_triples(:));
free_indices = all_indices(~ismember(psalm_names,excluded_psalms));

free_psalm_names_str = psalm_names(free_indices);
free_psalm_lengths = psalm_lengths(free_indices);

N_free = length(free_indices);

% Convert Psalm name notation from using trailing letters to a purely
% numeric notation (119a => 119.1). This allows us to calculate costs
% numerically.
free_psalm_names_num = zeros(N_free,1);
for i = 1:N_free
    free_psalm_names_num(i) = psalmNameStrToNum(free_psalm_names_str(i));
end


%% Candidate Generation
% Now that the data is in a workable state, we generate candidate triples.
combinations = nchoosek(1:N_free,3);
N_combinations = size(combinations,1);

% Next we filter the combinations based some of our constraints.
valid_indices = true(N_combinations,1);

for k = 1:N_combinations
    names = free_psalm_names_str(combinations(k,:));
    lengths = free_psalm_lengths(combinations(k,:));

    % Ensure that the triple is not too long or too short
    total_length = sum(lengths);
    if total_length > max_length || total_length < min_length
        valid_indices(k) = false;
        continue;
    end

    % Forbid multiple divisions of Psalm 119 from being in the same triple
    prefix_count = sum(startsWith(names, "119"));
    if prefix_count > 1
        valid_indices(k) = false;
        continue;
    end

    % Ensure that pairs are together
    for i = 1:numel(names)
        id = names(i);
        if isKey(pairs, id) && all(~ismember(names,pairs(id)))
            valid_indices(k) = false;
            break;
        end
    end

end

combinations = combinations(valid_indices,:);
N_combinations = size(combinations,1);


%% Cost Calculation
% Next, we need to calculate a "cost" for each triple. The lower the cost,
% the "better" the triple. This cost is a function of:
% a) how far away the length is from our desired length (the average)
% b) how close the Psalms are sequentially
cost = zeros(N_combinations,1);

for k = 1:N_combinations
    total_length = sum(free_psalm_lengths(combinations(k,:)));
    imbalance = abs(total_length - target_length);
    names = free_psalm_names_num(combinations(k,:));
    spread = max(names) - min(names);

    cost(k) = (w_val * imbalance) + (w_seq * spread);
end


%% Configure & Run Solver
% Now that we have built and filtered our list of qualified candidates, we
% can configure our solver variables based on our constraints and run the
% solver.

% Ensure that each Psalm is used exactly once
A = sparse(N_free, N_combinations);
for k = 1:N_combinations
    A(combinations(k,:), k) = 1;
end
Aeq = [A; ones(1,N_combinations)];

% Ensure that N/3 triples are selected
beq = [ones(N_free,1); N_free/3];

intcon = 1:N_combinations;
lb = zeros(N_combinations,1);
ub = ones(N_combinations,1);

result = intlinprog(cost, intcon, [], [], Aeq, beq, lb, ub);


%% Extract & Display Results
selected = find(result > 0.5);
solution = combinations(selected,:);
psalm_triples = free_indices(solution);
final_psalms = [];
final_lengths = [];

for i = 1:size(psalm_triples,1)
    t = psalm_triples(i,:);
    total_length = sum(psalm_lengths(t));
    final_psalms = [final_psalms; sort(psalm_names(t),'ascend')];
    final_lengths = [final_lengths; total_length];
end

for i = 1:size(forced_triples,1)
    t = forced_triples(i,:);
    total_length = ...
        psalm_lengths(psalm_names==t(1)) + ...
        psalm_lengths(psalm_names==t(2)) + ...
        psalm_lengths(psalm_names==t(3));
    final_psalms = [final_psalms; sort(t,'ascend')];
    final_lengths = [final_lengths; total_length];
end

T = table(final_psalms, final_lengths, VariableNames=["Psalms","Length"]);
T = sortrows(T,'Length','descend');


%% psalmNameStrToNum
% Helper function for converting the string name of a Psalm to a numeric
% one. E.g.: "119a" => 119.1
function n = psalmNameStrToNum(s)
    token = regexp(s,'(\d+)([a-z]?)','tokens');
    t = token{1};
    base = str2double(t{1});
    if isempty(t{2})
        offset = 0;
    else
        offset = double(t{2}) - double('a') + 1;
    end
    n = base + (offset * 0.1);
end

The original number of candidate triples is close to 600,000. Filtering reduces this to a little over 300,000 triples, which helps make the solver run much faster.

After running for about 10 minutes, the program produces the following optimal solution (optimal based on the constraints and weighting of goals).

Revisions and Final Result

These results are not bad, but the next step is to see if we can revise our constraints and goals to see if we can produce an even better result.

First, I removed all of the forced triples (except 78a/78b/78c). While I liked the ones I originally specified, they were all much shorter than the average triple length. Better to see if we can provide fewer constraints and just let the solver do its job.

The second improvement was to the cost function. Originally, the function applied a cost equal to the difference between the highest and lowest Psalm number times the weight w_seq. However, this does not adequetely capture our priorities. While we highly prioritize Psalms that are sequential, if there are non-sequential Psalms, we don’t really care how far apart they are. So I tried two revisions: in the first, we apply a penalty of 100 for each nonsequential pair of Psalms in a triplet (for a max possible penalty of 200), and then 0.1 cost for each step over 1 between each pair; in the second, we apply just the 100 penalty for nonsequential Psalms, and no 0.1 penalty per step. We compare the results between these approaches below.

The final improvement is to tighten the range of triple lengths. First to average+/-10, and then to average+/-7.5 and average+/-5. The latter two produced no results, so those ranges are too tight. But +/-10 did produce results, and is a nice improvement over +/-15.

Now that we have three different Psalm groupings, we can compare them. To do this, we can categorize the groupings into the following categories (each triple gets one category, with the following order of precedence).

For ordered triple [a b c]:

1. Perfectly Sequential: (b - a) <= 1 && (c - b) <= 1
2. Partial Division: round(a) == round(b) || round(b) == round(c)
3. Partially Sequential: (b - a) <= 1 || (c - b) <= 1
4. Assortment: everything else

We also show the standard deviation of all the lengths.

Of the four categories of triple, #4 “Assortment” is the least desirable. So we can immediately notice that our two revisions are vast improvements over the first attempt.

Next, between #2 “Partial Division” and #3 “Partially Sequential”, I do not think I can say that one is better than the other. Our two latter attempts only differ by one.

Also notable is that our second two attempts saw a big increase in category #1 “Perfectly Sequential”, from 14 to 31, over a 2x increase. This is by far our biggest improvement. The cost is the increase in Standard Deviation. Which I think is a cost worth paying.

Since our two revised attempts show almost identical results, including equivalent standard deviations for their lengths, it’s basically a wash, and it comes down to preference. Since, upon further reflection, I do not think the additional penalty of 0.1*(c-a) really provides much theoretical benefit (I do not necessarily see value in preferring slightly closer non-sequential Psalms, all else being equal), I have selected the last attempt as our final grouping. Here it is:

However, by “final”, I do not mean we are locked into this grouping. It is likely that we will want to make some small adjustments during the process of assigning groups to hours. But more on that later.

And for completeness, here is our final program:

N = 168;
psalm_names = ...
    ["89b", "119a", "22", "78a", "102", "118", "35", "73", "44", "18b", "136", ...
     "18a", "106b", "31", "119b", "119c", "119d", "119e", "119f", "119g", "78b", ...
     "55", "71", "139", "106a", "105b", "104a", "34", "38", "50", "74", "94", "105a", ...
     "107a", "37a", "25", "33", "103", "107b", "9", "49", "51", "77", "135", "145", ...
     "109a", "66", "72", "80", "147", "89a", "69b", "68a", "83", "88", "116", "37b", ...
     "69a", "7", "10", "40", "45", "59", "115", "132", "78c", "68b", "81", "86", "90", ...
     "91", "92", "17", "19", "48", "144", "21", "27", "39", "41", "56", ...
     "60", "65", "85", "108", "140", "148", "5", "30", "36", "62", "76", "79", "84", ...
     "96", "104b", "2", "26", "42", "46", "57", "58", "63", "97", "143", "109b", "6", ...
     "16", "29", "32", "52", "64", "75", "95", "8", "20", "24", "47", "111", "112", ...
     "141", "146", "3", "4", "12", "28", "54", "61", "98", "99", "113", "122", "137", ...
     "149", "67", "82", "101", "114", "121", "124", "129", "130", "138", "142", "11", ...
     "14", "53", "87", "110", "120", "1", "13", "23", "70", "126", "128", "150", ...
     "15", "43", "93", "100", "125", "127", "123", "131", "133", "134", "117"];
psalm_lengths = ...
    [34, 32, 32, 31, 29, 29, 28, 28, 27, 26, 26, 25, 25, 25, 24, 24, 24, 24, 24, 24, 24, 24, ...
     24, 24, 23, 23, 23, 23, 23, 23, 23, 23, 22, 22, 22, 22, 22, 22, 21, 21, 21, 21, 21, 21, ...
     21, 20, 20, 20, 20, 20, 19, 19, 19, 19, 19, 19, 18, 18, 18, 18, 18, 18, 18, 18, 18, 17, ...
     17, 17, 17, 17, 16, 16, 15, 15, 15, 15, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, ...
     13, 13, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 11, 11, 11, 11, ...
     11, 11, 11, 11, 11, 10, 10, 10, 10, 10, 10, 10, 10, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 8, ...
     8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 4, 3, 3, 3, 2];

% These parameters can be tweaked to produce different results.
w_val           = 1.0; % weight of objective that triple length should be close to average
non_seq_penalty = 100; % penalty for triples having non-sequential psalms
target_length   = mean(psalm_lengths) * 3; % target length of a triple of Psalms
max_length      = target_length + 10;      % maximum length of a triple of Psalms
min_length      = target_length - 10;      % minimum length of a triple of Psalms

% Force Psalms split into pairs to go together, codified as a key-value
% pair. Allows us to ensure that if we see <key> in a triple, that the
% triple also contains <val>.
pairs = dictionary( ...
    ["89a", "89b", "18a", "18b", "106a", "106b", "105a", "105b", "107a", "107b", ...
     "37a", "37b", "69a", "69b", "68a", "68b", "104a", "104b", "109a", "109b"], ...
    ["89b", "89a", "18b", "18a", "106b", "106a", "105b", "105a", "107b", "107a", ...
     "37b", "37a", "69b", "69a", "68b", "68a", "104b", "104a", "109b", "109a"]);

% Ensure that these triples of Psalms are grouped together.
forced_triples = [
                   "78a", "78b", "78c"
                 ];


%% Preprocessing
% The first step is to perform some simple transformations on the data so
% that it can more easily be processed.

% Exclude all "forced" Psalms from lists
all_indices     = 1 : N;
excluded_psalms = unique(forced_triples(:));
free_indices    = all_indices(~ismember(psalm_names, excluded_psalms));

free_psalm_names_str = psalm_names(free_indices);
free_psalm_lengths   = psalm_lengths(free_indices);

N_free = length(free_indices);

% Convert Psalm name notation from using trailing letters to a purely
% numeric notation (119a => 119.1). This allows us to calculate costs
% numerically.
free_psalm_names_num = zeros(N_free, 1);
for i = 1 : N_free
    free_psalm_names_num(i) = psalmNameStrToNum(free_psalm_names_str(i));
end


%% Candidate Generation
% Now that the data is in a workable state, we generate candidate triples.
combinations   = nchoosek(1 : N_free, 3);
N_combinations = size(combinations, 1);

% Next we filter the combinations based some of our constraints.
valid_indices = true(N_combinations, 1);

for k = 1 : N_combinations
    names   = free_psalm_names_str(combinations(k, :));
    lengths = free_psalm_lengths(combinations(k, :));

    % Ensure that the triple is not too long or too short
    total_length = sum(lengths);
    if total_length > max_length || total_length < min_length
        valid_indices(k) = false;
        continue;
    end

    % Forbid multiple divisions of Psalm 119 from being in the same triple
    prefix_count = sum(startsWith(names, "119"));
    if prefix_count > 1
        valid_indices(k) = false;
        continue;
    end

    % Ensure that pairs are together
    for i = 1 : numel(names)
        id = names(i);
        if isKey(pairs, id) && ~any(ismember(names, pairs(id)))
            valid_indices(k) = false;
            break;
        end
    end

end

combinations   = combinations(valid_indices, :);
N_combinations = size(combinations, 1);


%% Cost Calculation
% Next, we need to calculate a "cost" for each triple. The lower the cost,
% the "better" the triple. This cost is a function of:
% a) how far away the length is from our desired length (the average)
% b) how close the Psalms are sequentially
cost = zeros(N_combinations, 1);

for k = 1 : N_combinations
    total_length = sum(free_psalm_lengths(combinations(k, :)));
    imbalance    = abs(total_length - target_length);

    % Review the distance between the 1st/2nd and 2nd/3rd Psalm, sequentially.
    % Our priority is to give a big penalty to non-sequential Psalms.
    penalty = 0;
    names   = sort(free_psalm_names_num(combinations(k, :)), 'ascend');
    for i = 2 : 3
        difference = names(i) - names(i - 1);
        if difference > 1
            penalty = penalty + non_seq_penalty;
        end
    end
    
    cost(k) = (w_val * imbalance) + penalty;
end


%% Configure & Run Solver
% Now that we have built and filtered our list of qualified candidates, we
% can configure our solver variables based on our constraints and run the
% solver.

% Ensure that each Psalm is used exactly once
A = sparse(N_free, N_combinations);
for k = 1 : N_combinations
    A(combinations(k, :), k) = 1;
end
Aeq = [A; ones(1, N_combinations)];

% Ensure that N/3 triples are selected
beq = [ones(N_free, 1); N_free / 3];

intcon = 1 : N_combinations;
lb     = zeros(N_combinations, 1);
ub     = ones(N_combinations, 1);

result = intlinprog(cost, intcon, [], [], Aeq, beq, lb, ub);


%% Extract & Display Results
selected      = find(result > 0.5);
solution      = combinations(selected, :);
psalm_triples = free_indices(solution);
final_psalms  = [];
final_lengths = [];

for i = 1 : size(psalm_triples, 1)
    t            = psalm_triples(i, :);
    total_length = sum(psalm_lengths(t));
    final_psalms  = [final_psalms; sort(psalm_names(t), 'ascend')];
    final_lengths = [final_lengths; total_length];
end

for i = 1 : size(forced_triples, 1)
    t = forced_triples(i, :);
    total_length = ...
        psalm_lengths(psalm_names == t(1)) + ...
        psalm_lengths(psalm_names == t(2)) + ...
        psalm_lengths(psalm_names == t(3));
    final_psalms  = [final_psalms; sort(t, 'ascend')];
    final_lengths = [final_lengths; total_length];
end

T = table(final_psalms, final_lengths, VariableNames = ["Psalms", "Length"]);
T = sortrows(T, 'Length', 'descend');


%% psalmNameStrToNum
% Helper function for converting the string name of a Psalm to a numeric
% one. E.g.: "119a" => 119.1
function n = psalmNameStrToNum(s)
    token = regexp(s, '(\d+)([a-z]?)', 'tokens');
    t     = token{1};
    base  = str2double(t{1});
    if isempty(t{2})
        offset = 0;
    else
        offset = double(t{2}) - double('a') + 1;
    end
    n = base + (offset * 0.1);
end

Filling the Calendar

Comparing Psalm Placement in Existing Psalters

Before determining how we ought to place our 56 triples, we should look to see how our previously-reviewed Psalters place each Psalm. In addition, I have personally reviewed each Psalm to see if I have a particular inclination for how it should be placed. This is entirely personal preference. For example, for me, Psalm 51 (very penitential) is a Friday Psalm, as is Psalm 22 (prefiguring the Passion of our Lord). I have captured all this information in the following table:

Analyzing the Results

Next Steps and Future Work

Notes and References

• Schemas retrieved from http://www.gregorianbooks.com/gregorian/www/www.kellerbook.com/SCHEMA~1.HTM. These schemas here use the Greek numbering, which why it was preserved as the default on this page. The schemas also had incomplete data about which instances of Psalms were repetitions and which were divisions. An effort was made to add this data to the LotH table, but not to the others.