sunsets of a set 2

[ignacio_b]

Hello everyone,

I need to generate all the possible combinations of a set (in my case my set has 32 elements) and I need all possible subsets of 16 elements.
So, anyone have any idea how to do this?

There are subroutines and I have tried, however, documentation is not clear enough. So, if someone knows how to do this either by using a subroutine or by writing a code I would appreciate it.

Thanks and regards,

Ignacio

 

[xwb]

Have you tried the brute force method with 16 nested loops?

 

[ignacio_b]

Hello,

I have not tried yet, I am a bit new in FORTRAN... so, if you could give me an example it would help me a lot.

Thanks and kind regards,

 

[xwb]

Take a simple example 5 numbers 3 combinations. So that will be ₅C₃ where n = 5 and r = 3. Let us work it out manually. Always start with something small that you can work through manually. The combinations will be

1 2 3
1 2 4
1 2 5
1 3 4
1 3 5
1 4 5
2 3 4
2 3 5
2 4 5
3 4 5

A program will do exactly the same thing. Look at the first column. Notice how the first 6 start with 1, the next 3 start with 2 and the last one starts with 3. So your fist loop goes from 1 to (n - r + 1). Try it out

do i1 = 1, n - r + 1
    write(*,*) i1
end do

Now look at the second column. It goes from 2 to 4. i.e. from i1 + 1 to n - r + 2. Try it out

do i1 = 1, n - r + 1
    do i2 = i1 + 1, n - r + 2.
        write(*,*) i1, i2
    end do
end do

Now look at the third column. It goes from 3 to 5 i.e. from i2 + 1 to n - r + 3

So now you have your 3 loops.

do i1 = 1, n - r + 1
    do i2 = i1 + 1, n - r + 2
        do i3 = i2 + 1, n - r + 3
            write(*, *) i1, i2, i3
        end do
    end do
end do

That just prints the indices. OK now what happens if your numbers are not 1, 2, 3, 4, 5 but 3, 21, 16, 2, 5

integer, parameter:: n = 5, r = 3
integer:: num(n), comb(r)
! first stick them into an array
num = (/ 3, 21, 16, 2, 5 /)
! now pick out the indices
do i1 = 1, n - r + 1
    comb(1) = num(i1)
    do i2 = i1 + 1, n - r + 2
        comb(2) = num(i2)
        do i3 = i2 + 1, n - r + 3
            comb(3) = num(i3)
            write(*, *) (comb(ii), ii = 1, r, 1)
        end do
    end do
end do

Notice where the 1, 2 and 3 appear. There is a pattern there - I'll leave you to figure it out. Don't write the whole program in one go. Do a little bit at a time and then try it out. Once you are happy with it, move to the next stage.

 

[ignacio_b]

Hello,

First of all thank you so much for your reply and time.

I am going to check it and try to understand the meaning of these loops.
I know how to use loops but this is the very first time that I see this notation in loops.

Thanks again, I really appreciate it.

 

[ignacio_b]

Dear xwb's,

First of all, the code works perfectly and thanks for that. I just have some doubts due to my lack of experience in FORTRAN.

1.The 'i1' that you used in the loop, it is not specified in the variables section, and as far as I know it should be declared as an integer. It's clear that the code works but I don't get why you did not declared it.

2. Is there a way to storage every output (subset) into a separate file?

3. Could you give me a brief explanation of this: (comb(ii), ii = 1, r, 1). I am completely lost in this statement (of course, I know that you are printing out the results, but the reason of had used ii=1,r,1 is completely out of my understanding right now).

Sincerely,

Ignacio

 

[xwb]

1) By default, variables beginning with letters I to N are integer. They really ought to be declared but I was more concerned about going through the motions of developing a program than the syntax.
2) Yes, you need to look at the open, write and close statements. Try

https://www.tutorialspoint.com/fortran/fortran_file_input_output.htm

3) That is an implied do loop which prints elements comb(1) to comb(r). If you just look up implied do loop, there are lots of hits.

 

[mikrom]

I found nice iterative algorithm for generating combinations in lexicographic order on

https://www.baeldung.com/java-combinations-algorithm.

For example, for n = 5 and r = 3 it generates combinations starting from

1  2  3

in lexicographic order until

3  4  5

How it works: We set the initial combination = [1, 2, 3] and increase it by
combination(3) = combination(3) + 1 until combination(3) = 5, i.e. we generate:

1  2  3
1  2  4
1  2  5

now combination(3) = 5 so we change to combination(2) and increase it
combination(2) = combination(2) + 1 and set
combination(3) = combination(2) + 1 i.e.

1  3  4

and then combination(3) = combination(3) + 1 until combination(3) = 5, i.e.:

1  3  5

now combination(3) = 5 so we change to combination(2) and increase it
combination(2) = combination(2) + 1

1  4  5

now combination(3) = 5, combination(2) = 4 so we change to combination(1) and increase it by
combination(1) = combination(1) + 1 and for i = 2, 3 we set
combination(i) = combination(i - 1) + 1

2  3  4

then combination(3) = combination(3) + 1

2  3  5

now combination(3) = 5 again so we increase combination(2) = combination(2) + 1 and get

2  4  5

now combination(3) = 5, combination(2) = 4 so we change to combination(1) and increase it by
combination(1) = combination(1) + 1 and get finally

3  4  5

Here is the program:
combinations.f95

[COLOR=#800080]program[/color] combinations_test
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:, :), [COLOR=#2e8b57][b]allocatable[/b][/color] :: combinations
  [COLOR=#2e8b57][b]integer[/b][/color] :: n, r, nr_comb

  n [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]5[/color]
  r [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]3[/color]
  [COLOR=#0000ff]!n = 32[/color]
  [COLOR=#0000ff]!r = 16[/color]
  [COLOR=#008b8b]call[/color] combination_generator(n, r, [COLOR=#ff00ff].true.[/color])
[COLOR=#800080]end program[/color] combinations_test

[COLOR=#800080]subroutine[/color] combination_generator(n, r, print_combinations)
  [COLOR=#0000ff]! see:[/color]
  [COLOR=#0000ff]! [URL unfurl="true"]https://www.baeldung.com/java-combinations-algorithm[/URL][/color]
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: n
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: r
[COLOR=#2e8b57][b]  logical[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: print_combinations
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](r) :: combination
  [COLOR=#2e8b57][b]integer[/b][/color] :: t, i, k
  [COLOR=#6a5acd]10[/color] [COLOR=#a52a2a][b]format[/b][/color]([COLOR=#ff00ff]'Generated Combinations C('[/color],i2,[COLOR=#ff00ff]','[/color],i2,[COLOR=#ff00ff]'):'[/color])
  [COLOR=#6a5acd]20[/color] [COLOR=#a52a2a][b]format[/b][/color]([COLOR=#ff00ff]'C('[/color],i2,[COLOR=#ff00ff]','[/color],i2,[COLOR=#ff00ff]') ='[/color])

  [COLOR=#a52a2a][b]do[/b][/color] i[COLOR=#a52a2a][b]=[/b][/color][COLOR=#ff00ff]1[/color], r
    combination(i) [COLOR=#a52a2a][b]=[/b][/color] i
  [COLOR=#a52a2a][b]end do[/b][/color]

  [COLOR=#0000ff]! write header[/color]
  [COLOR=#a52a2a][b]if[/b][/color] (print_combinations) [COLOR=#a52a2a][b]then[/b][/color]
    [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#ff00ff]10[/color]) n, r
  [COLOR=#a52a2a][b]end if[/b][/color]    
  k [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  [COLOR=#a52a2a][b]do[/b][/color] [COLOR=#a52a2a][b]while[/b][/color](combination(r) [COLOR=#a52a2a][b]<=[/b][/color] n)
    k [COLOR=#a52a2a][b]=[/b][/color] k [COLOR=#a52a2a][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#a52a2a][b]if[/b][/color] (print_combinations) [COLOR=#a52a2a][b]then[/b][/color]
      [COLOR=#0000ff]! write one combination[/color]
      [COLOR=#a52a2a][b]do[/b][/color] i[COLOR=#a52a2a][b]=[/b][/color][COLOR=#ff00ff]1[/color], r
        [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#ff00ff]'(i3)'[/color],[COLOR=#a52a2a][b]advance[/b][/color][COLOR=#a52a2a][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) combination(i)
      [COLOR=#a52a2a][b]end do[/b][/color]
      [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color])
    [COLOR=#a52a2a][b]end if[/b][/color]

    [COLOR=#0000ff]! generate next combination in lexicographic order[/color]
    t [COLOR=#a52a2a][b]=[/b][/color] r
    [COLOR=#a52a2a][b]do[/b][/color] [COLOR=#a52a2a][b]while[/b][/color] ((t [COLOR=#a52a2a][b].ne.[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#a52a2a][b].and.[/b][/color] (combination(t) [COLOR=#a52a2a][b].eq.[/b][/color] n [COLOR=#a52a2a][b]-[/b][/color] r [COLOR=#a52a2a][b]+[/b][/color] t))
      t [COLOR=#a52a2a][b]=[/b][/color] t [COLOR=#a52a2a][b]-[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#a52a2a][b]end do[/b][/color]
    combination(t) [COLOR=#a52a2a][b]=[/b][/color] combination(t) [COLOR=#a52a2a][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#a52a2a][b]do[/b][/color] i [COLOR=#a52a2a][b]=[/b][/color] t [COLOR=#a52a2a][b]+[/b][/color] [COLOR=#ff00ff]1[/color], r
      combination(i) [COLOR=#a52a2a][b]=[/b][/color] combination(i [COLOR=#a52a2a][b]-[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#a52a2a][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#a52a2a][b]end do[/b][/color]
  [COLOR=#a52a2a][b]end do[/b][/color]
  [COLOR=#0000ff]! write nCr:     [/color]
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#ff00ff]20[/color], [COLOR=#a52a2a][b]advance[/b][/color][COLOR=#a52a2a][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) n, r
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) k
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color])
[COLOR=#800080]end subroutine[/color] combination_generator

You could get the resulting combination into the file using redirection:

$ ./combinations > combinations_out.txt

For bigger n, r use first rather

  [COLOR=#008b8b]call[/color] combination_generator(n, r, [COLOR=#ff00ff].false.[/color])

so the subroutine only computes the binomial coefficient C(n,r) and will not write the combinations, because writing large files will take very much time.
For example for n = 32, r = 16 it computes C(32,16) = 601080390. Creating the file with all combinations took very long time and the file has 28 GB. I cannot open it in editor, only display it's head and tail:

$ head combinations_32_16_out.txt 
Generated Combinations C(32,16):
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 17
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 18
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 19
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 20
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 21
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 22
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 23
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 24
$ tail combinations_32_16_out.txt 
 16 17 18 19 20 21 22 24 25 26 27 28 29 30 31 32
 16 17 18 19 20 21 23 24 25 26 27 28 29 30 31 32
 16 17 18 19 20 22 23 24 25 26 27 28 29 30 31 32
 16 17 18 19 21 22 23 24 25 26 27 28 29 30 31 32
 16 17 18 20 21 22 23 24 25 26 27 28 29 30 31 32
 16 17 19 20 21 22 23 24 25 26 27 28 29 30 31 32
 16 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
C(32,16) =   601080390

 

[ignacio_b]

Hi Mikrom,

Thanks for the detailed explanation, I appreciate it.

Regarding the file size, Yes it is huge. I need this code because I am writing a Monte Carlo Simulated Annealing (MC-SA) and I need to explore the arrangements of a particular crystal and each arrangement of 16 numbers [between 1-32] represent the indices of a set of 3D coordinates. Basically after generate an arrangement I need to pass it to another particular array, then explore/study it (By using MC-SA) and then do the same with the next array, so, technically speaking I don't need to store them in a file.

I've got one last question and it is about MPI: Is it possible to use MPI in a code like this?
As far as I know I would make the code run faster... Right?

I am sort of new in this "FORTRAN world", so I know how to implement some tasks but definitely some task surpass by much my skills... So, I am trying learn and understand as much as I can. So, If you know any book/manual or YouTube channel that you can recommend, it would be really helpful.

Once again, thank you so much.

Kind regards,

ignacio