DL_target_rate_calculator

LICENSE

Copyright 2020 Johannes Choo

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

rate.py

from typing import List, Tuple
from fractions import Fraction

# The compute_avg function in this script computes the average number of
# target 5*'s per pull on a banner, given
# 1. the banner's base 5* percent,
# 2. the pityless % of pulling a desired 5* unit, and
# 3. the number of pity-boosting single pulls you make before tenfolds
#
# This average count includes dupe targets from the same tenfold pull:
# that is, if you get two or more of the same target units in the same
# tenfold, the dupe units you get contribute to the average count, hence
# the average number of unique targets units you pull will be slightly
# lower.
#
# This average also assumes that your last pull (single/tenfold) on a
# banner includes a 5*. If you regularly let a banner rotate
# (hence resetting pity), your average is going to be slightly lower
# then the output average. Otherwise, your average will
# probabilistically converge to the average that this program computes.

# -----

# returns if the next pull should be a tenfold after num_pulled pulls
def is_tenfold(max_num_singles: int, num_pulled: int) -> bool:
  return num_pulled >= max_num_singles

# return the current_rate (for the next pull)
# of pulling a 5* unit given the number num_pulled of previous pulls,
# and the base (pitiless) rate of pulling a 5* unit of the banner.
# Historically, the pity rate always increases by 0.5% per 10 pulls,
# so we hardcode it here.
def current_rate(start_rate: Fraction, num_pulled: int) -> Fraction:
  return start_rate + Fraction(5, 1000) * (num_pulled // 10)

# given the current_rate of pulling a 5*, returns if the next pull
# contains a guaranteed five-star.
# Historically, this always happens when the current 5* rate is 9%,
# so we hardcode it in here.
def has_guaranteed(curr_rate: Fraction) -> bool:
  return curr_rate >= Fraction(9, 100)

# probability that next pull breaks pity if it is a single
def prob_broken_single(curr_rate: Fraction) -> Fraction:
  if has_guaranteed(curr_rate):
    return Fraction(1)
  else:
    return curr_rate

# average number of target units pulled in next pulls
# after num_pulled pulls if next pulls are a pity-breaking tenfold
#
# for pulls without a guaranteed pull,
# this follows a
# binomial distribution with zero-5* outcomes truncated away
# equivalent to conditional distribution given that outcome is not a
# zero-5* tenfold
# binomial mean is simply n * p = 10 * curr_rate
# probability of zero 5*'s pulled is
# binom(n, 0) * p**0 * (1-p)**(n-0) = (1-curr_rate) ** 10
# desired conditional expectation =
# (Sum(outcome in outcomes_that_are_not_outcome_without_5star) num_5star_in(outcome) * P(outcome)) / P(outcomes_that_are_not_outcome_without_5star)
# Since we have num_5star_in(outcome_without_5star) = 0, this is equal to
# (Sum(outcome in all_outcomes) num_5star_in(outcome) * P(outcome)) / P(outcomes_that_are_not_outcome_without_5star)
# = n * p / (1 - P(outcome_without_5star))
# = 10 * curr_rate / (1 - (1 - curr_rate) ** 10)
# then multiply this with target (probability that each 5* is a target unit) to obtain average number of target units pulled in this tenfold
# given that this tenfold contains a 5*.
#
# for pulls with a guaranteed pull, the probability of a 5* is 1 + 9 * curr_rate
def average_target_broken_tenfold(curr_rate: Fraction, target: Fraction) -> Fraction:
  if has_guaranteed(curr_rate):
    return target * (1 + 9 * curr_rate)
  else:
    return target * (10 * curr_rate) / (1 - (1 - curr_rate) ** 10)

# probability that next pulls break your pity if next pulls are
# a tenfold
def prob_broken_tenfold(curr_rate: Fraction) -> Fraction:
  if has_guaranteed(curr_rate):
    return Fraction(1)
  else:
    return 1 - (1 - curr_rate) ** 10

# Consider a round of pulls already until pity is broken. The round
# is already in progress, with num_pulled pulls already made, and only
# mass of all rounds get this far without being pity broken.
# gen_table(max_num_singles, start_rate, target, num_pulled, mass)
# given
# 1. max_num_singles is the number of single pulls before starting tenfold pulls
# 2. start_rate is the probability of pulling a 5* in a single with no pity
# 3. target is the probability that a 5* pulled in a single is a target unit
#    (this is independent of pity)
# 4. num_pulled is the number of pulls made so far
# 5. mass is the probability that a round has not terminated after
#    num_pulled pulls
# prints a list of (np, p, m) triples, where
# 1. np is the pulls made in this round before being pity broken,
# 2. p is the probability that a round is pity broken after exactly
#    np pulls
# 3. m is the average number of featured units pulled in this round
def gen_table(max_num_singles, start_rate, target, num_pulled: int, mass: Fraction) -> List[Tuple[int, Fraction, Fraction]]:
  if mass == 0:
    return []
  else:
    curr_rate = current_rate(start_rate, num_pulled)
    if is_tenfold(max_num_singles, num_pulled):
      next_num = num_pulled + 10
      p = mass * prob_broken_tenfold(curr_rate)
      m = average_target_broken_tenfold(curr_rate, target)
    else:
      next_num = num_pulled + 1
      p = mass * prob_broken_single(curr_rate)
      m = target
    l = gen_table(max_num_singles, start_rate, target, next_num, mass - p)
    l.append((next_num, p, m))
    return l

# 1. start_rate_percent is the percent probability of pulling a 5* on
# the banner without any
# pity. e.g. historically, non-gala had 4%, so put down 4, while
# gala had 6%, so put down 6
# 2. target_percent is the total percent probability of pulling
# your target unit (must be
# all 5* units). Note that listed percents in-game are usually the
# actual percent probabilty rounded down to several significant figures.
# If you want exact precision for a percent that cannot be represented
# exactly as a float, it's easy to modify this function to use Fraction
# instead.
def compute_avg(max_num_singles: int, start_rate_percent: float, target_percent: float):
  start_rate = Fraction(start_rate_percent) / 100
  target = Fraction(target_percent) / 100 / start_rate
  table = gen_table(max_num_singles, start_rate, target, 0, Fraction(1))
  average_targets_pulled = 0
  average_pulls_made = 0
  for num, p, m in table:
    average_targets_pulled += p * m
    average_pulls_made += p * num
  return average_targets_pulled / average_pulls_made

def generate_common_cases():
  maxes = [[0, 0.], [0, 0.], [0, 0.], [0, 0.], [0, 0.]]
  for i in range(102):
    ps = [
      float(compute_avg(i, 6, 0.5) * 100),
      float(compute_avg(i, 4, 0.5) * 100),
      float(compute_avg(i, 4, 1) * 100),
      float(compute_avg(i, 4, 1.5) * 100),
      float(compute_avg(i, 4, 2) * 100)
    ]
    for j in range(len(ps)):
      if ps[j] > maxes[j][1]:
        maxes[j] = [i, ps[j]]
    print("{}: 6%, .5%: {:.6}%; 4%, .5%: {:.6}%; 4%, 1%: {:.6}%; 4%, 1.5%: {:.6}%; 4%, 2%: {:.6}%".format(
      i,
      ps[0],
      ps[1],
      ps[2],
      ps[3],
      ps[4]
    ))
  print("6%, .5%: {}|{:.6}%; 4%, .5%: {}|{:.6}%; 4%, 1%: {}|{:.6}%; 4%, 1.5%: {}|{:.6}%; 4%, 2%: {}|{:.6}%".format(
    maxes[0][0], maxes[0][1],
    maxes[1][0], maxes[1][1],
    maxes[2][0], maxes[2][1],
    maxes[3][0], maxes[3][1],
    maxes[4][0], maxes[4][1],
  ))

# unused function. prints the statistics among the round cases we've
# considered so far as we compute the average by going through round
# cases.
def print_table(max_num_singles: int, start_rate_percent: float, target_percent: float):
  start_rate = Fraction(start_rate_percent) / 100
  target = Fraction(target_percent) / 100 / start_rate
  table = gen_table(max_num_singles, start_rate, target, 0, Fraction(1))
  average_targets_pulled = 0
  average_pulls_made = 0
  for num, p, m in reversed(table):
    average_targets_pulled += p * m
    average_pulls_made += p * num
    print("{:>3}, p: {:<8.5}, this_targs: {:<8.5}, cumtarg: {:<8.5}, cumpulls: {}, avg: {}".format(
      num,
      float(p), # fraction of rounds that broke in previous single / tenfold
      float(m), # average number of target units pulled in last pull
      float(average_targets_pulled), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_pulls_made), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_targets_pulled / average_pulls_made) # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      ),
    )

generate_common_cases()

rate_dupeless.py

from typing import Dict, Iterable, List, Tuple
from fractions import Fraction

from sympy.ntheory import multinomial_coefficients

# The compute_avg function in this script computes the average number of
# unique target 5*'s per pull on a banner, given
# a distribution (generated by get_distr) of the relevant outcomes
#
# This average count counts dupe targets from the same tenfold pull only
# once:
# that is, if you get two or more of the same target units in the same
# tenfold, the dupe units you get do not contribute to the average
# count, hence the average number of targets units you pull will
# be slightly higher.
#
# This average also assumes that your last pull (single/tenfold) on a
# banner includes a 5*. If you regularly let a banner rotate
# (hence resetting pity), your average is going to be slightly lower
# then the output average. Otherwise, your average will
# probabilistically converge to the average that this program computes.

# -----

def ninefold_coefficients_closure():
  memo: Dict[int, Dict[tuple, int]] = {}
  def ninefold_coefficients(m: int) -> Dict[tuple, int]:
    if m in memo:
      return dict(memo[m])
    else:
      memo[m] = multinomial_coefficients(m, 9)
      return dict(memo[m])
  return ninefold_coefficients
ninefold_coefficients = ninefold_coefficients_closure()

def tenfold_coefficients_closure():
  memo: Dict[int, Dict[tuple, int]] = {}
  def tenfold_coefficients(m) -> Dict[tuple, int]:
    if m in memo:
      return dict(memo[m])
    else:
      memo[m] = multinomial_coefficients(m, 10)
      return dict(memo[m])
  return tenfold_coefficients
tenfold_coefficients = tenfold_coefficients_closure()

# Construct a probability distribution that we recognize from the
# following information
# 1. start_rate_up_adv_permille is the permille (i.e. times thousand)
#    probability of pulling any rate up 5*
#    adventurer with a single pull on the banner without any
#    pity. Typically (5 * number of rate up adventurers)
# 2. start_rate_up_drag_permille is as start_rate_up_adv_permille, but
#    for rate up dragons instead. Typically (8 * number of rate up dragons)
# 3. start_rest_adv_permille is as start_rate_up_adv_permille, but for
#    non-rate-up adventurers. Typically (20 - start_rate_up_adv_permille
#    for non-gala, 30 - start_rate_up_adv_permille for gala)
# 4. start_rest_drag_permille is as start_rest_adv_permille, but for
#    non-rate-up dragons. Typically (20 - start_rate_up_drag_permille
#    for non-gala, 30 - start_rate_up_drag_permille for gala)
# 5. num_rate_up_adv is the number of rate up adventurers on the banner
# 6. num_rate_up_drag is the number of rate up dragons on the banner
# 7. num_rest_adv is the number of non-rate up adventurers on the banner
# 8. num_rest_drag is the number of non-rate up dragons on the banner
# 9. num_target_rate_up_adv is the number of targeted rate up
#    adventurers on the banner
# 10. num_target_rate_up_drag is the number of targeted rate up
#     dragons on the banner
# 11. num_rest_rate_up_adv is the number of targeted non rate up
#    adventurers on the banner
# 12. num_rest_rate_up_drag is the number of targeted non rate up
#    dragons on the banner
def get_distr(
  start_rate_up_adv_permille: float,
  start_rate_up_drag_permille: float,
  start_rest_adv_permille: float,
  start_rest_drag_permille: float,
  num_rate_up_adv: int,
  num_rate_up_drag: int,
  num_rest_adv: int,
  num_rest_drag: int,
  num_target_rate_up_adv: int,
  num_target_rate_up_drag: int,
  num_target_rest_adv: int,
  num_target_rest_drag: int
) -> List[Fraction]:
  starts: Tuple[Fraction, ...] = tuple(map(lambda i: Fraction(i) / 1000, [
    start_rate_up_adv_permille, start_rate_up_drag_permille, start_rest_adv_permille, start_rest_drag_permille
  ]))
  nums: Tuple[int, ...] = (num_rate_up_adv, num_rate_up_drag, num_rest_adv, num_rest_drag)
  num_targets: Tuple[int, ...] = (num_target_rate_up_adv, num_target_rate_up_drag, num_target_rest_adv, num_target_rest_drag)
  start_pers: Tuple[Fraction, ...] = tuple(starts[i] / nums[i] if starts[i] > 0 else 0 for i in range(4))
  start_non_5star: Fraction = Fraction(1) - sum(starts)
  start_non_target: Fraction = Fraction(sum(starts[i] - start_pers[i] * num_targets[i] for i in range(4)))
  distr: List[Fraction] = [start_non_5star, start_non_target]
  labels = ["non_5star", "non_target"]
  for i in range(num_target_rate_up_adv):
    distr.append(start_pers[0])
    labels.append("rate_up_adv_" + str(i))
  for i in range(num_target_rate_up_drag):
    distr.append(start_pers[1])
    labels.append("rate_up_drag_" + str(i))
  for i in range(num_target_rest_adv):
    distr.append(start_pers[2])
    labels.append("rest_adv_" + str(i))
  for i in range(num_target_rest_drag):
    distr.append(start_pers[3])
    labels.append("rest_drag_" + str(i))
  return distr

# returns if the next pull should be a tenfold after num_pulled pulls
def is_tenfold(max_num_singles: int, num_pulled: int) -> bool:
  return num_pulled >= max_num_singles

# Adds pity to the current distribution
# Historically, the pity rate always increases by 0.5% per 10 pulls,
# so we hardcode it here.
def add_pity(distr: List[Fraction]) -> None:
  non_5star = distr[0]
  five_star = 1 - non_5star
  for i in range(len(distr)):
    distr[i] = distr[i] + Fraction(5, 1000) * distr[i] / five_star
  distr[0] = non_5star - Fraction(5, 1000)

# given the current_rate of pulling a 5*, returns if the next pull
# contains a guaranteed five-star.
# Historically, this always happens when the current 5* rate is 9%,
# so we hardcode it in here.
def has_guaranteed(distr: List[Fraction]) -> bool:
  return 1 - distr[0] >= Fraction(9, 100)

# probability that next pull breaks pity if it is a single
def prob_broken_single(distr: List[Fraction]) -> Fraction:
  if has_guaranteed(distr):
    return Fraction(1)
  else:
    return 1 - distr[0]

# probability that next pull pulls a target if it is a single
# that pulls a 5*
def value_single(distr: List[Fraction]) -> Fraction:
  return (1 - distr[0] - distr[1]) / (1 - distr[0])

# average number of unique target units pulled in next tenfold pulls
def value_tenfold(distr: List[Fraction]) -> Fraction:
  if has_guaranteed(distr):
    cond_distr = distr[1:]
    no_5star_single = distr[0]
    five_star_single = 1 - distr[0]
    # scale cond_distr
    for i in range(len(cond_distr)):
      cond_distr[i] /= five_star_single
    # cond_distr now contains the distribution for the guaranteed 5* pull
    # we generate the multinomial distribution for the ninefold
    multi_distr = ninefold_coefficients(len(distr))
    # for each multinomial case
    for k in multi_distr:
      for i in range(len(k)):
        multi_distr[k] *= distr[i] ** k[i]
    # multi_distr now contains the multinomial distribution for the ninefold
    value = Fraction(0)
    for i in range(len(cond_distr)):
      for k in multi_distr:
        # sum all valuable (i.e. nonzeros on target units) nonzeros
        # k_value is the value of this kind of tenfold pull.
        # The slice excludes counting the non-5star pulls and the non-target-5star pulls
        # if there are any such.
        k_value = sum(1 for j in range(2, len(k)) if k[j] > 0 or i + 1 == j)
        # use below formula instead if we want to count dupes
        # k_value = sum(k[j] + (1 if i + 1 == j else 0) for j in range(2, len(k)) if k[j] > 0 or i + 1 == j)
        # we now scale the value by the probability of this case appearing in a guaranteed + ninefold
        # then add the scaled value to what eventually becomes the average value of a guaranteed + ninefold,
        # as the sum of all the different cases scaled by probability of appearance.
        value += multi_distr[k] * cond_distr[i] * k_value
    return value
  else:
    # easier to directly compute the value in multi_distr corresponding to key (10, 0, 0, ..., 0)
    no_5star = distr[0] ** 10
    five_star = 1 - no_5star
    multi_distr = tenfold_coefficients(len(distr))
    for k in multi_distr:
      # compute joint probability for particular multinomial case
      for i in range(len(k)):
        multi_distr[k] *= distr[i] ** k[i]
      # normalize against probability of pulling at least a 5star to get conditional probability
      multi_distr[k] /= five_star
    value = Fraction(0)
    # luckily, case we normalized over must contribute zero to value in below summation
    # (summation of value per multinomial case)
    for k in multi_distr:
      # k_value is the value of this kind of tenfold pull.
      # The slice excludes counting the non-5star pulls and the non-target-5star pulls
      # if there are any such.
      k_value = sum(1 for ki in k[2:] if ki > 0)
      # use below formula instead if we want to count dupes
      # k_value = sum(ki for ki in k[2:] if ki > 0)
      # add this value to total value with contribution modified by probability of occurrence
      # we now scale the value by the probability of this case appearing in a tenfold
      # then add the scaled value to what eventually becomes the average value of a tenfold,
      # as the sum of all the different cases scaled by probability of appearance.
      value += multi_distr[k] * k_value
    return value

# probability that next pulls break your pity if next pulls are
# a tenfold
def prob_broken_tenfold(distr: List[Fraction]) -> Fraction:
  if has_guaranteed(distr):
    return Fraction(1)
  else:
    return 1 - distr[0] ** 10

# gen_table(max_num_singles, distr)
# given
# 1. max_num_singles is the number of single pulls before starting tenfold pulls
# 2. distr is the distribution of a single pull (at the start)
# prints a list of (np, p, m) triples, where
# 1. np is the pulls made in this round before being pity broken,
# 2. p is the probability that a round is pity broken after exactly
#    np pulls
# 3. m is the average value of 5* units pulled in this round
# Note: distr is mutated by this function
def gen_table(max_num_singles, distr: List[Fraction]) -> List[Tuple[int, Fraction, Fraction]]:
  num_pulled = 0
  mass = Fraction(1)
  table: List[Tuple[int, Fraction, Fraction]] = []
  while mass != 0:
    if is_tenfold(max_num_singles, num_pulled):
      num_pulled += 10
      p = mass * prob_broken_tenfold(distr)
      table.append((num_pulled, p, value_tenfold(distr)))
      mass -= p
      add_pity(distr)
    else:
      num_pulled += 1
      p = mass * prob_broken_single(distr)
      table.append((num_pulled, p, value_single(distr)))
      mass -= p
      if num_pulled % 10 == 0:
        add_pity(distr)
  return table

# given a table generated by gen_table
# computes the average value per pull, averaging over all the cases of rounds
def compute_avg(table: List[Tuple[int, Fraction, Fraction]]) -> Fraction:
  average_targets_pulled = Fraction(0)
  average_pulls_made = Fraction(0)
  for num, p, m in table:
    average_targets_pulled += p * m
    average_pulls_made += p * num
  return average_targets_pulled / average_pulls_made

# compute and print the average value per pull across different numbers of single summons
# if targeting only the gala adventurer
"""
def print_gala_table():
  distr = get_distr(
    start_rate_up_adv_permille=5,
    start_rate_up_drag_permille=0,
    start_rest_adv_permille=25,
    start_rest_drag_permille=30,
    num_rate_up_adv=1,
    num_rate_up_drag=0,
    num_rest_adv=20, # unimportant if not targeting non-rate-up advs
    num_rest_drag=20, # unimportant if not targeting non-rate-up dragons
    num_target_rate_up_adv=1,
    num_target_rate_up_drag=0,
    num_target_rest_adv=0,
    num_target_rest_drag=0
  )
  table = gen_table(0, distr)
  average_targets_pulled = Fraction(0)
  average_pulls_made = Fraction(0)
  cum_p = Fraction(0)
  for num, p, m in table:
    average_targets_pulled += p * m
    average_pulls_made += p * num
    cum_p += p
    print("{:>3}, p: {:<8.5}, cum_p: {:<8.5}, this_targs: {:<8.5}, cumtarg: {:<8.5}, cumpulls: {}, avg: {}".format(
      num,
      float(p), # fraction of rounds that broke in previous single / tenfold
      float(cum_p), # fraction of rounds that broke in previous single / tenfold
      float(m), # average number of target units pulled in last pull
      float(average_targets_pulled), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_pulls_made), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_targets_pulled / average_pulls_made) # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      ),
    )
  print(float(compute_avg(table)))
"""

# compute and print the average value per pull across different numbers of single summons
# if targeting only the gala adventurer
def compute_gala():
  maxes = (0, 0.)
  for i in range(102):
    distr = get_distr(
      start_rate_up_adv_permille=5,
      start_rate_up_drag_permille=0,
      start_rest_adv_permille=25,
      start_rest_drag_permille=30,
      num_rate_up_adv=1,
      num_rate_up_drag=0,
      num_rest_adv=20, # wrong, but unimportant if not targeting non-rate-up advs
      num_rest_drag=20, # wrong, but unimportant if not targeting non-rate-up dragons
      num_target_rate_up_adv=1,
      num_target_rate_up_drag=0,
      num_target_rest_adv=0,
      num_target_rest_drag=0
    )
    table = gen_table(i, distr)
    curr_avg = float(compute_avg(table))
    print("{}: 6%, .5%: {:.6}%".format(
      i,
      curr_avg * 100
    ))
    if curr_avg > maxes[1]:
      maxes = (i, curr_avg)
  print("max: {}: 6%, .5%: {:.6}%".format(maxes[0], maxes[1] * 100))

# compute and print the average value per pull across different numbers of single summons
# if targeting only the rate-up adventurers on the FEH banner
def compute_feh():
  maxes = (0, 0.)
  for i in range(102):
    distr = get_distr(
      start_rate_up_adv_permille=15,
      start_rate_up_drag_permille=0,
      start_rest_adv_permille=10,
      start_rest_drag_permille=15,
      num_rate_up_adv=3,
      num_rate_up_drag=0,
      num_rest_adv=9+10+11+7+11,
      num_rest_drag=7+7+9+7+6,
      num_target_rate_up_adv=3,
      num_target_rate_up_drag=0,
      num_target_rest_adv=0,
      num_target_rest_drag=0
    )
    table = gen_table(i, distr)
    curr_avg = float(compute_avg(table))
    print("{}: FEH all rate-up: {:.6}%".format(
      i,
      curr_avg * 100
    ))
    if curr_avg > maxes[1]:
      maxes = (i, curr_avg)
  print("max: {}: FEH all rate-up: {:.6}%".format(maxes[0], maxes[1] * 100))

compute_feh()

rate_value.py

from typing import Dict, Iterable, List, Tuple
from fractions import Fraction

from sympy.ntheory import multinomial_coefficients

# The compute_avg function in this script computes the average number of
# unique target 5*'s per pull on a banner, given
# a distribution (generated by get_distr) of the relevant outcomes
#
# This average count counts dupe targets from the same tenfold pull only
# once:
# that is, if you get two or more of the same target units in the same
# tenfold, the dupe units you get do not contribute to the average
# count, hence the average number of targets units you pull will
# be slightly higher.
#
# This average also assumes that your last pull (single/tenfold) on a
# banner includes a 5*. If you regularly let a banner rotate
# (hence resetting pity), your average is going to be slightly lower
# then the output average. Otherwise, your average will
# probabilistically converge to the average that this program computes.

# -----

def ninefold_coefficients_closure():
  memo: Dict[int, Dict[tuple, int]] = {}
  def ninefold_coefficients(m: int) -> Dict[tuple, int]:
    if m in memo:
      return dict(memo[m])
    else:
      memo[m] = multinomial_coefficients(m, 9)
      return dict(memo[m])
  return ninefold_coefficients
ninefold_coefficients = ninefold_coefficients_closure()

def tenfold_coefficients_closure():
  memo: Dict[int, Dict[tuple, int]] = {}
  def tenfold_coefficients(m) -> Dict[tuple, int]:
    if m in memo:
      return dict(memo[m])
    else:
      memo[m] = multinomial_coefficients(m, 10)
      return dict(memo[m])
  return tenfold_coefficients
tenfold_coefficients = tenfold_coefficients_closure()

# Construct a probability distribution that we recognize from the
# following information
# 1. start_rate_up_adv_permille is the permille (i.e. times thousand)
#    probability of pulling any rate up 5*
#    adventurer with a single pull on the banner without any
#    pity. Typically (5 * number of rate up adventurers)
# 2. start_rate_up_drag_permille is as start_rate_up_adv_permille, but
#    for rate up dragons instead. Typically (8 * number of rate up dragons)
# 3. start_rest_adv_permille is as start_rate_up_adv_permille, but for
#    non-rate-up adventurers. Typically (20 - start_rate_up_adv_permille
#    for non-gala, 30 - start_rate_up_adv_permille for gala)
# 4. start_rest_drag_permille is as start_rest_adv_permille, but for
#    non-rate-up dragons. Typically (20 - start_rate_up_drag_permille
#    for non-gala, 30 - start_rate_up_drag_permille for gala)
# 5. num_rate_up_adv is the number of rate up adventurers on the banner
# 6. num_rate_up_drag is the number of rate up dragons on the banner
# 7. num_rest_adv is the number of non-rate up adventurers on the banner
# 8. num_rest_drag is the number of non-rate up dragons on the banner
# 9. num_target_rate_up_adv is the number of targeted rate up
#    adventurers on the banner
# 10. num_target_rate_up_drag is the number of targeted rate up
#     dragons on the banner
# 11. num_rest_rate_up_adv is the number of targeted non rate up
#    adventurers on the banner
# 12. num_rest_rate_up_drag is the number of targeted non rate up
#    dragons on the banner
def get_distr(
  start_rate_up_adv_permille: float,
  start_rate_up_drag_permille: float,
  start_rest_adv_permille: float,
  start_rest_drag_permille: float,
  num_rate_up_adv: int,
  num_rate_up_drag: int,
  num_rest_adv: int,
  num_rest_drag: int,
  target_rate_up_adv: List[List[Fraction]],
  target_rate_up_drag: List[List[Fraction]],
  target_rest_adv: List[List[Fraction]],
  target_rest_drag: List[List[Fraction]]
) -> Tuple[List[Fraction], List[List[Fraction]]]:
  starts: Tuple[Fraction, ...] = tuple(map(lambda i: Fraction(i) / 1000, [
    start_rate_up_adv_permille, start_rate_up_drag_permille, start_rest_adv_permille, start_rest_drag_permille
  ]))
  nums: Tuple[int, ...] = (num_rate_up_adv, num_rate_up_drag, num_rest_adv, num_rest_drag)
  num_target: Tuple[int, ...] = (len(target_rate_up_adv), len(target_rate_up_drag), len(target_rest_adv), len(target_rest_drag))
  start_per: Tuple[Fraction, ...] = tuple(starts[i] / nums[i] if starts[i] > 0 else 0 for i in range(4))
  start_non_5star: Fraction = Fraction(1) - sum(starts)
  start_non_target: Fraction = Fraction(sum(starts[i] - start_per[i] * num_target[i] for i in range(4)))
  distr: List[Fraction] = [start_non_5star, start_non_target]
  value: List[List[Fraction]] = [[], []]
  # initialize distr and values with ["non_5star", "non_target"]
  for unit_values in target_rate_up_adv:
    distr.append(start_per[0])
    value.append(unit_values)
  for unit_values in target_rate_up_drag:
    distr.append(start_per[1])
    value.append(unit_values)
  for unit_values in target_rest_adv:
    distr.append(start_per[2])
    value.append(unit_values)
  for unit_values in target_rest_drag:
    distr.append(start_per[3])
    value.append(unit_values)
  return distr, value

# returns if the next pull should be a tenfold after num_pulled pulls
def is_tenfold(max_num_singles: int, num_pulled: int) -> bool:
  return num_pulled >= max_num_singles

# Adds pity to the current distribution
# Historically, the pity rate always increases by 0.5% per 10 pulls,
# so we hardcode it here.
def add_pity(distr: List[Fraction]) -> None:
  non_5star = distr[0]
  five_star = 1 - non_5star
  for i in range(len(distr)):
    distr[i] = distr[i] + Fraction(5, 1000) * distr[i] / five_star
  distr[0] = non_5star - Fraction(5, 1000)

# given the current_rate of pulling a 5*, returns if the next pull
# contains a guaranteed five-star.
# Historically, this always happens when the current 5* rate is 9%,
# so we hardcode it in here.
def has_guaranteed(distr: List[Fraction]) -> bool:
  return 1 - distr[0] >= Fraction(9, 100)

# probability that next pull breaks pity if it is a single
def prob_broken_single(distr: List[Fraction]) -> Fraction:
  if has_guaranteed(distr):
    return Fraction(1)
  else:
    return 1 - distr[0]

# probability that next pull pulls a target if it is a single
# that pulls a 5*
def value_single(distr: List[Fraction], unit_value: List[List[Fraction]]) -> Fraction:
  value = Fraction(0)
  for i in range(len(distr)):
    if unit_value[i]:
      value += distr[i] * unit_value[i][0]
  return value / (1 - distr[0])

# average number of unique target units pulled in next tenfold pulls
def value_tenfold(distr: List[Fraction], unit_value: List[List[Fraction]]) -> Fraction:
  if has_guaranteed(distr):
    # we generate the multinomial distribution for the ninefold
    multi_distr = ninefold_coefficients(len(distr))
    # for each multinomial case
    for k in multi_distr:
      # for each kind
      for i in range(len(k)):
        # scale distribution case by the probability that k[i] units of the kind occurs in the case
        multi_distr[k] *= distr[i] ** k[i]
    # multi_distr now contains the multinomial distribution for the ninefold
    value = Fraction(0)
    for i in range(1, len(distr)):
      for k in multi_distr:
        # sum all valuable (i.e. nonzeros on target units) nonzeros
        # k_value is the value of this kind of tenfold pull.
        # The slice excludes counting the non-5star pulls and the non-target-5star pulls
        # if there are any such.
        k_value: Fraction = sum(
          # for each kind j, add the per-dupe (indexed by n) value to k_value
          (sum((unit_value[j][n] for n in range(min(k[j] + (1 if i == j else 0), len(unit_value[j])))), Fraction(0)) for j in range(len(k))),
          Fraction(0))
        # use below formula instead if we want to count dupes
        # k_value = sum(k[j] + (1 if i + 1 == j else 0) for j in range(2, len(k)) if k[j] > 0 or i + 1 == j)
        # we now scale the value by the probability of this case appearing in a guaranteed + ninefold
        # then add the scaled value to what eventually becomes the average value of a guaranteed + ninefold,
        # as the sum of all the different cases scaled by probability of appearance.
        value += multi_distr[k] * distr[i] * k_value
    return value / (1 - distr[0])
  else:
    multi_distr = tenfold_coefficients(len(distr))
    for k in multi_distr:
      # compute joint probability for particular multinomial case
      for i in range(len(k)):
        multi_distr[k] *= distr[i] ** k[i]
    value = Fraction(0)
    # luckily, case we normalized over must contribute zero to value in below summation
    # (summation of value per multinomial case)
    for k in multi_distr:
      # k_value is the value of this kind of tenfold pull.
      # The slice excludes counting the non-5star pulls and the non-target-5star pulls
      # if there are any such.
      k_value = sum(
        (sum((unit_value[j][n] for n in range(min(k[j], len(unit_value[j])))), Fraction(0)) for j in range(len(k))),
        Fraction(0))
      # use below formula instead if we want to count dupes
      # k_value = sum(ki for ki in k[2:] if ki > 0)
      # add this value to total value with contribution modified by probability of occurrence
      # we now scale the value by the probability of this case appearing in a tenfold
      # then add the scaled value to what eventually becomes the average value of a tenfold,
      # as the sum of all the different cases scaled by probability of appearance.
      value += multi_distr[k] * k_value
    # easier to directly compute the value in multi_distr corresponding to key (10, 0, 0, ..., 0)
    return value / (1 - distr[0] ** 10)

# probability that next pulls break your pity if next pulls are
# a tenfold
def prob_broken_tenfold(distr: List[Fraction]) -> Fraction:
  if has_guaranteed(distr):
    return Fraction(1)
  else:
    return 1 - distr[0] ** 10

# gen_table(max_num_singles, distr)
# given
# 1. max_num_singles is the number of single pulls before starting tenfold pulls
# 2. distr is the distribution of a single pull (at the start)
# prints a list of (np, p, m) triples, where
# 1. np is the pulls made in this round before being pity broken,
# 2. p is the probability that a round is pity broken after exactly
#    np pulls
# 3. m is the average value of 5* units pulled in this round
# Note: distr is mutated by this function
def gen_table(max_num_singles, distr: List[Fraction], value: List[List[Fraction]]) -> List[Tuple[int, Fraction, Fraction]]:
  num_pulled = 0
  mass = Fraction(1)
  table: List[Tuple[int, Fraction, Fraction]] = []
  while mass != 0:
    if is_tenfold(max_num_singles, num_pulled):
      num_pulled += 10
      p = mass * prob_broken_tenfold(distr)
      table.append((num_pulled, p, value_tenfold(distr, value)))
      mass -= p
      add_pity(distr)
    else:
      num_pulled += 1
      p = mass * prob_broken_single(distr)
      table.append((num_pulled, p, value_single(distr, value)))
      mass -= p
      if num_pulled % 10 == 0:
        add_pity(distr)
  return table

# given a table generated by gen_table
# computes the average value per pull, averaging over all the cases of rounds
def compute_avg(table: List[Tuple[int, Fraction, Fraction]]) -> Fraction:
  average_targets_pulled = Fraction(0)
  average_pulls_made = Fraction(0)
  for num, p, m in table:
    average_targets_pulled += p * m
    average_pulls_made += p * num
  return average_targets_pulled / average_pulls_made

# compute and print the average value per pull across different numbers of single summons
# if targeting only the gala adventurer
def print_gala_table():
  distr, value = get_distr(
    start_rate_up_adv_permille=15,
    start_rate_up_drag_permille=0,
    start_rest_adv_permille=10,
    start_rest_drag_permille=15,
    num_rate_up_adv=3,
    num_rate_up_drag=0,
    num_rest_adv=9+10+11+7+11,
    num_rest_drag=7+7+9+7+6,
    target_rate_up_adv=[[Fraction(1)], [Fraction(1)], [Fraction(1)]],
    target_rate_up_drag=[],
    target_rest_adv=[],
    target_rest_drag=[]
  )
  table = gen_table(20, distr, value)
  average_targets_pulled = Fraction(0)
  average_pulls_made = Fraction(0)
  cum_p = Fraction(0)
  for num, p, m in table:
    average_targets_pulled += p * m
    average_pulls_made += p * num
    cum_p += p
    print("{:>3}, p: {:<8.5}, cum_p: {:<8.5}, this_targs: {:<8.5}, cumtarg: {:<8.5}, cumpulls: {}, avg: {}".format(
      num,
      float(p), # fraction of rounds that broke in previous single / tenfold
      float(cum_p), # fraction of rounds that broke in previous single / tenfold
      float(m), # average number of target units pulled in last pull
      float(average_targets_pulled), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_pulls_made), # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      float(average_targets_pulled / average_pulls_made) # average number of target units pulled per pull among rounds that stop at this number of pulls or more
      ),
    )
  print(float(compute_avg(table)))

# compute and print the average value per pull across different numbers of single summons
# if targeting only the gala adventurer
def compute_gala():
  target_rate_up_adv: List[List[Fraction]] = []
  target_rate_up_drag: List[List[Fraction]] = [
    [Fraction(1, 1), Fraction(1, 8), Fraction(1, 8), Fraction(1, 8)] # Mars
  ]
  target_rest_adv: List[List[Fraction]] = []
  target_rest_drag: List[List[Fraction]] = []
  maxes = (0, 0.)
  for i in range(102):
    distr, value = get_distr(
      start_rate_up_adv_permille=5,
      start_rate_up_drag_permille=0,
      start_rest_adv_permille=25,
      start_rest_drag_permille=30,
      num_rate_up_adv=0,
      num_rate_up_drag=1,
      num_rest_adv=20, # wrong, but unimportant if not targeting non-rate-up advs
      num_rest_drag=20, # wrong, but unimportant if not targeting non-rate-up dragons
      target_rate_up_adv=target_rate_up_adv,
      target_rate_up_drag=target_rate_up_drag,
      target_rest_adv=target_rest_adv,
      target_rest_drag=target_rest_drag
    )
    table = gen_table(i, distr, value)
    curr_avg = float(compute_avg(table))
    print("{}: 6%, .5%: {:.6}%".format(
      i,
      curr_avg * 100
    ))
    if curr_avg > maxes[1]:
      maxes = (i, curr_avg)
  print("max: {}: 6%, .5%: {:.6}%".format(maxes[0], maxes[1] * 100))

print_gala_table()