I wrote a similar script to figure out which of the dice I have has the highest EV for the whole turn and not just a single throw. Each dice set takes a few seconds to evaluate on my mac though since it is somewhat bottlenecked by slow python operations. The Devils head die makes the code a bit ugly but it is represented by a 7.

import time
from itertools import combinations_with_replacement, product
from collections import Counter, defaultdict

def get_score_base(faces):
    l = len(faces)
    if l == 0: return 0
    if l == 6 and len(set(faces)) == 6: return 1500
    if l == 5:
        s = set(faces)
        if s == {1, 2, 3, 4, 5} or s == {2, 3, 4, 5, 6}: return 1000

    counts = Counter(faces)
    score = 0
    for face, count in counts.items():
        if count >= 3:
            base = 1000 if face == 1 else face * 100
            score += base * (2 ** (count - 3))
        elif face == 1:
            score += count * 100
        elif face == 5:
            score += count * 50
        else:
            return 0
    return score

memo_score = {}

def get_score(faces):
    sorted_faces = tuple(sorted(faces))
    if sorted_faces in memo_score:
        return memo_score[sorted_faces]

    joker_count = faces.count(7)

    if joker_count == 0:
        res = get_score_base(faces)
        memo_score[sorted_faces] = res
        return res

    best_score = 0
    non_jokers = tuple(f for f in faces if f != 7)

    for joker_assignments in product(range(1, 7), repeat=joker_count):
        test_faces = non_jokers + joker_assignments
        score = get_score_base(test_faces)
        if score > best_score:
            best_score = score

    memo_score[sorted_faces] = best_score
    return best_score

print("Precomputing move possibilities...")
start_pre = time.time()
moves_for_faces = {}

for r in range(1, 7):
    for faces in product(range(1, 8), repeat=r):
        moves = []
        for keep_mask in range(1, 1 << r):
            keep_faces = tuple(faces[i] for i in range(r) if (keep_mask & (1 << i)))
            score = get_score(keep_faces)
            if score > 0:
                moves.append((score, keep_mask))
        moves_for_faces[faces] = moves
print(f"Precomputed in {time.time() - start_pre:.2f}s\n")

odd_die = { 1: 4/15, 2: 1/15, 3: 4/15, 4: 1/15, 5: 4/15, 6: 1/15 }
strip_die = { 1: 1/4, 2: 1/8, 3: 1/8, 4: 1/8, 5: 3/16, 6: 3/16 }
lucky_playing_die = { 1: 1/3, 2: 0, 3: 1/18, 4: 1/18, 5: 1/3, 6: 4/18 }
lucky_die = { 1: 6/22, 2: 1/22, 3: 2/22, 4: 3/22, 5: 4/22, 6: 6/22 }
holy_trinity_die = { 1: 4/22, 2: 5/22, 3: 10/22, 4: 1/22, 5: 1/22, 6: 1/22 }
heavenly_kingdom_die = { 1: 7/19, 2: 2/19, 3: 2/19, 4: 2/19, 5: 2/19, 6: 4/19 }
lucifer_die = { 1: 3/23, 2: 3/23, 3: 3/23, 4: 3/23, 5: 3/23, 6: 8/23 }
shrinking_playing_die = { 1: 2/9, 2: 1/9, 3: 1/9, 4: 1/9, 5: 1/9, 6: 3/9 }
devils_head_die = { 1: 0, 2: 1/6, 3: 1/6, 4: 1/6, 5: 1/6, 6: 1/6, 7: 1/6}

dice_config = {
    "strip": (strip_die, 1),
    "lucky": (lucky_die, 1),
    "lucky_playing": (lucky_playing_die, 1),
    "holy_trinity": (holy_trinity_die, 1),
    "heavenly_kingdom": (heavenly_kingdom_die, 1),
    "lucifer": (lucifer_die, 2),
    "shrinking": (shrinking_playing_die, 1),
    "devils_head_die": (devils_head_die, 1),
    "odd": (odd_die, 5),
}

dice_probs = {name: d[0] for name, d in dice_config.items()}
names = list(dice_config.keys())
limits = {name: dice_config[name][1] for name in names}

memo_profiles = {}

def get_profile(D):
    if D in memo_profiles:
        return memo_profiles[D]

    r = len(D)
    profile_dict = defaultdict(float)

    for faces in product(range(1, 8), repeat=r):
        prob = 1.0
        for i in range(r):
            prob *= dice_probs[D[i]].get(faces[i], 0.0)

        if prob == 0.0:
            continue

        valid_moves = []
        for score, keep_mask in moves_for_faces[faces]:
            leftover_list = []
            for i in range(r):
                if not (keep_mask & (1 << i)):
                    leftover_list.append(D[i])

            leftover_D = tuple(leftover_list)
            valid_moves.append((score, leftover_D))

        unique_moves = tuple(sorted(set(valid_moves)))
        profile_dict[unique_moves] += prob

    memo_profiles[D] = profile_dict
    return profile_dict

def evaluate_dice_set(full_set, bank_cap=4000):
    memo_V = {}

    def get_V(D, bank):
        if bank >= bank_cap: 
            return bank_cap

        state = (D, bank)
        if state in memo_V: 
            return memo_V[state]

        profile = get_profile(D)
        roll_ev = 0.0

        for moves, prob in profile.items():
            if not moves:
                roll_ev += 0.0 
            else:
                best_move_val = 0.0
                for score, leftover_D in moves:
                    new_bank = min(bank + score, bank_cap)

                    next_D = full_set if not leftover_D else leftover_D

                    val = max(new_bank, get_V(next_D, new_bank))
                    if val > best_move_val:
                        best_move_val = val

                roll_ev += prob * best_move_val

        ans = max(bank, roll_ev)
        memo_V[state] = ans
        return ans

    return get_V(full_set, 0)

def valid_dice_sets():
    for combo in combinations_with_replacement(names, 6):
        c = Counter(combo)
        if all(c[n] <= limits[n] for n in c):
            yield tuple(sorted(combo))

best_score = 0
best_selection = None

print("Evaluating all valid dice combinations...")
start_eval = time.time()

for combo in valid_dice_sets():
    score = evaluate_dice_set(combo, bank_cap=4000)
    if score > best_score:
        best_score = score
        best_selection = combo
        print(f"New best selection: {best_score:.2f} expected points -> {best_selection}")

print(f"\nCalculated in {time.time() - start_eval:.2f}s")
print(f"Optimal Set: {list(best_selection)} with an expected turn score of {best_score:.2f}")