Digits Dilemma

By: Jonathan Carroll

Re-posted from: https://jcarroll.com.au/2024/06/15/digits-dilemma/

Another day, another short riddle to be solved with several programming
languages! This one is nice because solving it doesn’t need a lot of code, but
it uses some interesting aspects of evaluation.

I saw this post
which isn’t new (it’s from 2022) that poses a nice problem to solve:

With the numbers 123456789, make them add up to 100.
They must stay in the same order but you can use addition, subtraction,
multiplication, division, brackets etc. All numbers must be used exactly once.

and demonstrates a solution in Haskell

import Control.Monad (forM)
import Language.Haskell.Interpreter

expressions :: [String]
expressions =
  let ops = [ '+', '-', '/', '*' ]
   in [ [ '1', a, '2', b, '3', c, '4', d, '5', e, '6', f, '7', g, '8', h, '9' ]
        | a <- ops, b <- ops
        , c <- ops, d <- ops
        , e <- ops, f <- ops
        , g <- ops, h <- ops
        ]

result = runInterpreter $ do
  setImports ["Prelude"]
  exprs <- forM expressions evaluate
  pure $ filter (\(_, a) -> a == "100") $ fromRight [] exprs
  where

I’m trying to learn Haskell this year, so it was a great opportunity try to
follow along. I’m still working my way towards being able to run short snippets
– the ghci tool for interactive use has a bit of a learning curve and doesn’t
immediately let me use the imports (or I’m doing something wrong) so while I
think I can follow the steps as presented, I can’t yet dig into them as
interactively as I’d like.

The general idea, though, is to use a comprehension to expand all combinations of
the allowed operators (+, -, /, and *) between the values 1 to 9. I’m
somewhat familiar with comprehensions and played with them in my post
Pythagorean Triples with Comprehensions
in several languages, including Haskell.

I wanted to see how I might go about this problem in R, and I knew I’d need to make
some adjustments because R does not have comprehensions.

One way to get all the combinations of operators between the values is to use
expand.grid() which generates all combinations of its inputs

expand.grid(1:3, letters[1:3])
##   Var1 Var2
## 1    1    a
## 2    2    a
## 3    3    a
## 4    1    b
## 5    2    b
## 6    3    b
## 7    1    c
## 8    2    c
## 9    3    c

Defining the operators as strings, I can generate a data.frame of the values and
all combinations of operators between them

ops <- c("*", "+", "-", "/")
combos <- expand.grid(1, ops, 2, ops, 3, ops, 4, ops, 5, ops, 6, ops, 7, ops, 8, ops, 9)
head(combos)
##   Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8 Var9 Var10 Var11 Var12 Var13 Var14
## 1    1    *    2    *    3    *    4    *    5     *     6     *     7     *
## 2    1    +    2    *    3    *    4    *    5     *     6     *     7     *
## 3    1    -    2    *    3    *    4    *    5     *     6     *     7     *
## 4    1    /    2    *    3    *    4    *    5     *     6     *     7     *
## 5    1    *    2    +    3    *    4    *    5     *     6     *     7     *
## 6    1    +    2    +    3    *    4    *    5     *     6     *     7     *
##   Var15 Var16 Var17
## 1     8     *     9
## 2     8     *     9
## 3     8     *     9
## 4     8     *     9
## 5     8     *     9
## 6     8     *     9

This generates a lot of combinations – with 4 possible operators in 8 possible
positions there are \(4^8\) = 65,536 combinations.

Pasting these numbers and operators together into expressions

exprs <- apply(combos, 1, \(x) paste0(x, collapse = ""))
head(exprs)
## [1] "1*2*3*4*5*6*7*8*9" "1+2*3*4*5*6*7*8*9" "1-2*3*4*5*6*7*8*9"
## [4] "1/2*3*4*5*6*7*8*9" "1*2+3*4*5*6*7*8*9" "1+2+3*4*5*6*7*8*9"

I get something I can evaluate as if I typed 1*2*3 into a console. I can get
the results of evaluating those with

results <- sapply(exprs, \(x) eval(parse(text = x)))
head(results)
## 1*2*3*4*5*6*7*8*9 1+2*3*4*5*6*7*8*9 1-2*3*4*5*6*7*8*9 1/2*3*4*5*6*7*8*9 
##            362880            362881           -362879             90720 
## 1*2+3*4*5*6*7*8*9 1+2+3*4*5*6*7*8*9 
##            181442            181443

Now, I just need to see which of those produces a value of 100. Because sapply
produced a vector with the expression itself as the name, I can extract the names
of the results which are equal to 100

answers <- names(which(results == 100))
answers
##  [1] "1*2*3-4*5+6*7+8*9" "1+2+3-4*5+6*7+8*9" "1+2-3*4-5+6*7+8*9"
##  [4] "1-2*3-4-5+6*7+8*9" "1+2-3*4+5*6+7+8*9" "1-2*3-4+5*6+7+8*9"
##  [7] "1-2*3+4*5+6+7+8*9" "1*2*3+4+5+6+7+8*9" "1+2+3+4+5+6+7+8*9"
## [10] "1+2*3+4*5-6+7+8*9" "1+2*3*4*5/6+7+8*9" "1*2*3*4+5+6-7+8*9"
## [13] "1*2*3*4+5+6+7*8+9" "1-2+3*4*5-6+7*8-9" "1-2+3*4*5+6*7+8-9"

Any of those can easily be verified manually

1*2*3*4+5+6+7*8+9
## [1] 100

One thing I noticed here was that I have one more result than the Haskell post
produces

length(answers)
## [1] 15

One of the answers stands out in that it contains a division, and sure enough this
is the one that doesn’t appear in the Haskell post. I’m not quite sure why – I
think the operator precedence is the same between R and Haskell, at least in terms
of these expressions

3 / 2 + 1
## [1] 2.5
ghci> 3 / 2 + 1
2.5

But since I haven’t yet been able to actually run that Haskell code myself, I can’t
verify those solutions.

My R solution to this puzzle is then

ops <- c("*", "+", "-", "/")
combos <- expand.grid(1, ops, 2, ops, 3, ops, 4, ops, 5, ops, 6, ops, 7, ops, 8, ops, 9)
exprs <- apply(combos, 1, \(x) paste0(x, collapse = ""))
results <- sapply(exprs, \(x) eval(parse(text = x)))
names(which(results == 100))

I did want to try a language which does have comprehensions – let’s try Julia!

Setting up the comprehension makes for a bit of a long line, but comes out okay

ops = ['+', '-', '*', '/']
## 4-element Vector{Char}:
##  '+': ASCII/Unicode U+002B (category Sm: Symbol, math)
##  '-': ASCII/Unicode U+002D (category Pd: Punctuation, dash)
##  '*': ASCII/Unicode U+002A (category Po: Punctuation, other)
##  '/': ASCII/Unicode U+002F (category Po: Punctuation, other)
exprs = ['1' * a * '2' * b * '3' * c * '4' * d * '5' * e * '6' * f * '7' * g * '8' * h * '9' 
         for a in ops, b in ops, c in ops, 
             d in ops, e in ops, f in ops, 
             g in ops, h in ops];
first(exprs, 10)
## 10-element Vector{String}:
##  "1+2+3+4+5+6+7+8+9"
##  "1-2+3+4+5+6+7+8+9"
##  "1*2+3+4+5+6+7+8+9"
##  "1/2+3+4+5+6+7+8+9"
##  "1+2-3+4+5+6+7+8+9"
##  "1-2-3+4+5+6+7+8+9"
##  "1*2-3+4+5+6+7+8+9"
##  "1/2-3+4+5+6+7+8+9"
##  "1+2*3+4+5+6+7+8+9"
##  "1-2*3+4+5+6+7+8+9"

That produces an array with many dimensions

size(exprs)
## (4, 4, 4, 4, 4, 4, 4, 4)

so it needs to be flattened into a vector with vec(). From there, it’s similar
to the R approach and I can use a eval(Meta.parse()) pattern, keeping in mind
that one can ‘broadcast’ scalar operations to vector operations with the dot
(.) operator

results = eval.(Meta.parse.(vec(exprs)));
first(results, 10)
## 10-element Vector{Real}:
##  45
##  41
##  44
##  42.5
##  39
##  35
##  38
##  36.5
##  46
##  34

Finding the values equal to 100 is similar to the R approach

exprs[findall(results .== 100)]
## 15-element Vector{String}:
##  "1*2*3*4+5+6+7*8+9"
##  "1-2+3*4*5+6*7+8-9"
##  "1-2+3*4*5-6+7*8-9"
##  "1+2+3+4+5+6+7+8*9"
##  "1*2*3+4+5+6+7+8*9"
##  "1-2*3+4*5+6+7+8*9"
##  "1+2*3+4*5-6+7+8*9"
##  "1-2*3-4+5*6+7+8*9"
##  "1+2-3*4+5*6+7+8*9"
##  "1+2*3*4*5/6+7+8*9"
##  "1*2*3*4+5+6-7+8*9"
##  "1-2*3-4-5+6*7+8*9"
##  "1+2-3*4-5+6*7+8*9"
##  "1+2+3-4*5+6*7+8*9"
##  "1*2*3-4*5+6*7+8*9"

and again we see the 15 answers including the one with a division, confirming the
R result.

This was a fun exploration – I don’t think I would want to try to solve it without
code, but the code solutions were a great opportunity to use a few different
languages.

I suspect there are a few different ways of solving this apart from this brute-force
expanding of every combination, perhaps with a solver or something. If you have one, I’d
love to see it. I thought @coolbutuseless had done something like this but the closest
I could find was this post
which is slightly different.

If you have comments, suggestions, or improvements, as always, feel free to use
the comment section below, or hit me up on
Mastodon.

devtools::session_info()
 
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.3.3 (2024-02-29)
##  os       Pop!_OS 22.04 LTS
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_AU.UTF-8
##  ctype    en_AU.UTF-8
##  tz       Australia/Adelaide
##  date     2024-07-06
##  pandoc   3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date (UTC) lib source
##  blogdown      1.18    2023-06-19 [1] CRAN (R 4.3.2)
##  bookdown      0.36    2023-10-16 [1] CRAN (R 4.3.2)
##  bslib         0.6.1   2023-11-28 [3] CRAN (R 4.3.2)
##  cachem        1.0.8   2023-05-01 [3] CRAN (R 4.3.0)
##  callr         3.7.3   2022-11-02 [3] CRAN (R 4.2.2)
##  cli           3.6.1   2023-03-23 [1] CRAN (R 4.3.3)
##  crayon        1.5.2   2022-09-29 [3] CRAN (R 4.2.1)
##  devtools      2.4.5   2022-10-11 [1] CRAN (R 4.3.2)
##  digest        0.6.34  2024-01-11 [3] CRAN (R 4.3.2)
##  ellipsis      0.3.2   2021-04-29 [3] CRAN (R 4.1.1)
##  evaluate      0.23    2023-11-01 [3] CRAN (R 4.3.2)
##  fastmap       1.1.1   2023-02-24 [3] CRAN (R 4.2.2)
##  fs            1.6.3   2023-07-20 [3] CRAN (R 4.3.1)
##  glue          1.7.0   2024-01-09 [1] CRAN (R 4.3.3)
##  htmltools     0.5.7   2023-11-03 [3] CRAN (R 4.3.2)
##  htmlwidgets   1.6.2   2023-03-17 [1] CRAN (R 4.3.2)
##  httpuv        1.6.12  2023-10-23 [1] CRAN (R 4.3.2)
##  icecream      0.2.1   2023-09-27 [1] CRAN (R 4.3.2)
##  jquerylib     0.1.4   2021-04-26 [3] CRAN (R 4.1.2)
##  jsonlite      1.8.8   2023-12-04 [3] CRAN (R 4.3.2)
##  JuliaCall     0.17.5  2022-09-08 [1] CRAN (R 4.3.3)
##  knitr         1.45    2023-10-30 [3] CRAN (R 4.3.2)
##  later         1.3.1   2023-05-02 [1] CRAN (R 4.3.2)
##  lifecycle     1.0.4   2023-11-07 [1] CRAN (R 4.3.3)
##  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.3.3)
##  memoise       2.0.1   2021-11-26 [3] CRAN (R 4.2.0)
##  mime          0.12    2021-09-28 [3] CRAN (R 4.2.0)
##  miniUI        0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
##  pkgbuild      1.4.2   2023-06-26 [1] CRAN (R 4.3.2)
##  pkgload       1.3.3   2023-09-22 [1] CRAN (R 4.3.2)
##  prettyunits   1.2.0   2023-09-24 [3] CRAN (R 4.3.1)
##  processx      3.8.3   2023-12-10 [3] CRAN (R 4.3.2)
##  profvis       0.3.8   2023-05-02 [1] CRAN (R 4.3.2)
##  promises      1.2.1   2023-08-10 [1] CRAN (R 4.3.2)
##  ps            1.7.6   2024-01-18 [3] CRAN (R 4.3.2)
##  purrr         1.0.2   2023-08-10 [3] CRAN (R 4.3.1)
##  R6            2.5.1   2021-08-19 [1] CRAN (R 4.3.3)
##  Rcpp          1.0.11  2023-07-06 [1] CRAN (R 4.3.2)
##  remotes       2.4.2.1 2023-07-18 [1] CRAN (R 4.3.2)
##  rlang         1.1.4   2024-06-04 [1] CRAN (R 4.3.3)
##  rmarkdown     2.25    2023-09-18 [3] CRAN (R 4.3.1)
##  rstudioapi    0.15.0  2023-07-07 [3] CRAN (R 4.3.1)
##  sass          0.4.8   2023-12-06 [3] CRAN (R 4.3.2)
##  sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.3.2)
##  shiny         1.7.5.1 2023-10-14 [1] CRAN (R 4.3.2)
##  stringi       1.8.3   2023-12-11 [3] CRAN (R 4.3.2)
##  stringr       1.5.1   2023-11-14 [3] CRAN (R 4.3.2)
##  urlchecker    1.0.1   2021-11-30 [1] CRAN (R 4.3.2)
##  usethis       2.2.2   2023-07-06 [1] CRAN (R 4.3.2)
##  vctrs         0.6.5   2023-12-01 [1] CRAN (R 4.3.3)
##  xfun          0.41    2023-11-01 [3] CRAN (R 4.3.2)
##  xtable        1.8-4   2019-04-21 [1] CRAN (R 4.3.2)
##  yaml          2.3.8   2023-12-11 [3] CRAN (R 4.3.2)
## 
##  [1] /home/jono/R/x86_64-pc-linux-gnu-library/4.3
##  [2] /usr/local/lib/R/site-library
##  [3] /usr/lib/R/site-library
##  [4] /usr/lib/R/library
## 
## ──────────────────────────────────────────────────────────────────────────────

Coin-tossing game: a numerical approach

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2024/06/14/probability3.html

Introduction

Today I decided to follow up on my last post solving a coin-tossing game.
This time, instead of simulation I want to use numerical approach
(and so probably a bit harder).

The post was written under Julia 1.10.1 and Graphs.jl 1.11.0.

The problem

Let me describe the setting of a game first (it is an extension of this post).

Assume Alice and Bob toss a fair coin. In each toss head (h) or tail (t) can show up with equal probability.

Alice and Bob choose some sequence of h and t they are waiting for. We assume that the chosen sequences have the same length and are different.
For example Alice could choose htht and Bob tthh.

The winner of the game is the person who saw their sequence first.

The question we ask if for a fixed sequence length n we can get cycles, that is, for example, that sequence s1 beats s2, s2 beats s3, and s3 beats s1.

To answer this question we will represent the game as a Markov process.

Step 1: non-terminating Markov chain

First we create a transition matrix of a Markov chain tracking current n element sequence in the game we consider.
Here is the code:

function markov(size::Integer)
    idx2states = vec(join.(Iterators.product([['h', 't'] for _ in 1:size]...)))
    states2idx = Dict(idx2states .=> eachindex(idx2states))
    P = zeros(2^size, 2^size)
    for state in idx2states
        for next in ("h", "t")
            nextstate = chop(state, head=1, tail=0) * next
            P[states2idx[state], states2idx[nextstate]] = 0.5
        end
    end
    return P, idx2states, states2idx
end

What we do in it is as follows:

  1. idx2states vector keeps track of all h and t sequences that have length n (i.e. it is a mapping from state number to state signature).
  2. states2idx is an inverse mapping – from state signature to state number.
  3. P is transition matrix of our chain. Note that from the sequence ab... (where all elements are h or t) we go to sequence b...h or b...t with equal probability.

Step 2: terminating Markov chain

We now need to create a function that is aware of Alice’s and Bob’s chosen sequences and make them terminating. We want to compute the probabilities of ending up in Alice’s and Bob’s state.
Here is the code:

function game(P, states2idx, alice, bob)
    P_game = copy(P)
    alice_idx, bob_idx = states2idx[alice], states2idx[bob]
    P_game[alice_idx, :] .= 0.0
    P_game[alice_idx, alice_idx] = 1.0
    P_game[bob_idx, :] .= 0.0
    P_game[bob_idx, bob_idx] = 1.0
    n = length(states2idx)
    terminal = fill(1 / n, 1, n) * P_game^(2^30)
    return terminal[states2idx[alice]], terminal[states2idx[bob]]
end

Note that we first update the P_game matrix to make alice_idx and bob_idx states terminating. Then, since I was lazy, we assume we make 2^30 steps of the process (fortunately in Julia it is fast).
Observe that initially all states are equally probably, so terminal matrix keeps information about long term probabilities of staying in all possible states.
We extract the probabilities of Alice’s and Bob’s states and return them.

Step 3: looking for cycles

We are now ready for a final move. We can consider all possible preferred sequences of Alice and Bob and create a graph that keeps track of which sequences beat other sequences:

using Graphs

function analyze_game(size::Integer, details::Bool=true)
    P, idx2states, states2idx = markov(size)
    g = SimpleDiGraph(length(states2idx))
    details && println("\nWinners:")
    for alice in idx2states, bob in idx2states
        alice > bob || continue
        alice_win, bob_win = game(P, states2idx, alice, bob)
        if alice_win > 0.51
            winner = "alice"
            add_edge!(g, states2idx[alice], states2idx[bob])
        elseif bob_win > 0.51
            winner = "bob"
            add_edge!(g, states2idx[bob], states2idx[alice])
        else
            winner = "tie (or close :))"
        end
        details && println(alice, " vs ", bob, ": ", winner)
    end
    cycles = simplecycles(g)
    if !isempty(cycles)
        min_len = minimum(length, cycles)
        filter!(x -> length(x) == min_len, cycles)
    end
    println("\nCycles:")
    for cycle in cycles
        println(idx2states[cycle])
    end
end

Note that I used 0.51 threshold for detection of dominance of one state over the other. We could do it better, but in practice for small n it is enough and working this way is simpler numerically.
What this threshold means is that we want to be “sure” that one player beats the other.
In our code we do two things:

  • optionally print the information which state beats which state;
  • print information about cycles found in beating patterns (we keep only cycles of the shortest length).

Let us check the code. Start with sequences of length 2:

julia> analyze_game(2)

Winners:
th vs hh: alice
th vs ht: tie (or close :))
ht vs hh: tie (or close :))
tt vs hh: tie (or close :))
tt vs th: tie (or close :))
tt vs ht: bob

Cycles:

We see that only th beats hh and ht beats tt (this is a symmetric case). We did not find any cycles.

Let us check 3:

julia> analyze_game(3)

Winners:
thh vs hhh: alice
thh vs hth: tie (or close :))
thh vs hht: alice
thh vs htt: tie (or close :))
hth vs hhh: alice
hth vs hht: bob
tth vs hhh: alice
tth vs thh: alice
tth vs hth: alice
tth vs hht: tie (or close :))
tth vs tht: alice
tth vs htt: bob
hht vs hhh: tie (or close :))
tht vs hhh: alice
tht vs thh: tie (or close :))
tht vs hth: tie (or close :))
tht vs hht: bob
tht vs htt: tie (or close :))
htt vs hhh: alice
htt vs hth: tie (or close :))
htt vs hht: bob
ttt vs hhh: tie (or close :))
ttt vs thh: bob
ttt vs hth: bob
ttt vs tth: tie (or close :))
ttt vs hht: bob
ttt vs tht: bob
ttt vs htt: bob

Cycles:
["thh", "hht", "htt", "tth"]

We now have the cycle. The shortest cycle has length 4 and it is unique. Let us see what happens for patterns of length 4 (I suppress printing the details as there are too many of them):

julia> analyze_game(4, false)

Cycles:
["thhh", "hhth", "hthh"]
["thhh", "hhtt", "ttth"]
["hhth", "hthh", "thht"]
["hhth", "thtt", "tthh"]
["hthh", "hhtt", "thth"]
["hthh", "hhtt", "ttht"]
["hthh", "hhtt", "ttth"]
["thht", "hhtt", "ttth"]
["htht", "thtt", "tthh"]
["thtt", "tthh", "hhht"]
["thtt", "htth", "ttht"]
["thtt", "httt", "ttht"]
["tthh", "hhht", "htth"]
["tthh", "hhht", "httt"]

In this case we have many cycles that are even shorter as they have length three.

Conclusions

The conclusion is that the game is slightly surprising. We can have cycles of dominance between sequences. I hope you liked this example. Happy summer!

JuliaHub Air: Secure, High-Performance Computing for Government Innovation

By: Jasmine Chokshi

Re-posted from: https://info.juliahub.com/blog/high-performance-computing-for-government-innovation

Secure, high-performance computing with JuliaHub Air for government innovation.

Government agencies rely heavily on High-Performance Computing (HPC) for a wide range of critical tasks, including scientific research, simulations, data analysis, and decision support. However, deploying and managing HPC solutions in secure, air-gapped environments presents unique challenges related to data protection, compliance, and restricted network access.