API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Mar 18th, 2016
381
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 42.32 KB | None | 0 0
raw download clone embed print report
  1. !  RUBIK.f90
  2. !
  3. !  FUNCTIONS:
  4. !   RUBIK      - Entry point of console application.
  5. !
  6. !
  7.  
  8. !****************************************************************************
  9. !
  10. !  PROGRAM: RUBIK
  11. !
  12. !
  13. !****************************************************************************
  14.  
  15. module external_fcts
  16. implicit none
  17.  
  18. integer  :: m, n, o
  19.  
  20. contains
  21.  
  22.     subroutine ToString(sequence_numbers, numbersToString)
  23.  
  24.         integer, dimension(:), intent(in)   :: sequence_numbers
  25.         integer                             :: a, taille_sequence
  26.         character, dimension(:),intent(out) :: numbersToString
  27.         taille_sequence = size(sequence_numbers)
  28.        
  29.  
  30.         do a = 1, taille_sequence
  31.             select case(mod(sequence_numbers(a),6))
  32.                 case(0)
  33.                     numbersToString(3*a-2) = 'F'
  34.                 case(1)
  35.                     numbersToString(3*a-2) = 'U'
  36.                 case(2)
  37.                     numbersToString(3*a-2) = 'R'
  38.                 case(3)
  39.                     numbersToString(3*a-2) = 'B'
  40.                 case(4)
  41.                     numbersToString(3*a-2) = 'D'
  42.                 case(5)
  43.                     numbersToString(3*a-2) = 'L'
  44.             end select
  45.            
  46.             numbersToString(3*a-1) = ' '
  47.             if (sequence_numbers(a)> 5) numbersToString(3*a-1) = '2'
  48.             if (sequence_numbers(a)>11) numbersToString(3*a-1) = "'"
  49.             numbersToString(3*a) = ' '
  50.         end do
  51.  
  52.     end subroutine ToString
  53.  
  54.     subroutine ToNumbers(sequence_string, StringToNumbers)
  55.         character, dimension(2), intent(in)     :: sequence_string
  56.         integer, intent(out)                    :: StringToNumbers
  57.         integer                                 :: a, b, taille_string
  58.  
  59.         do a=1, 2
  60.             select case(sequence_string(a))
  61.                 case('F')
  62.                     StringToNumbers = 0
  63.                 case('U')
  64.                     StringToNumbers = 1
  65.                 case('R')
  66.                     StringToNumbers = 2
  67.                 case('B')
  68.                     StringToNumbers = 3
  69.                 case('D')
  70.                     StringToNumbers = 4
  71.                 case('L')
  72.                     StringToNumbers = 5
  73.                 case("2")
  74.                     StringToNumbers = StringToNumbers+6
  75.                 case("'")
  76.                     StringToNumbers = StringToNumbers+12
  77.                 !case default
  78.                 !   print*, 'Incorrect caracter. Exit and run the program again'
  79.  
  80.             end select
  81.         end do
  82.         !print*, StringToNumbers
  83.  
  84.     end subroutine ToNumbers
  85.  
  86.     function init(cube)
  87.         integer, dimension(0:53), intent(in)  :: cube
  88.         integer, dimension(0:53)              :: init
  89.         integer                               :: indice = 0, temp = 0
  90.        
  91.         do m = 0, 5
  92.             do n = 0, 8
  93.             init(indice) = temp
  94.             indice = indice+1
  95.             end do
  96.             temp = temp+1
  97.         end do
  98.  
  99.     return
  100.     end function init
  101.  
  102.  
  103.     function do_move(cube_state, move)
  104.         integer, dimension(0:53), intent(in) :: cube_state
  105.         integer, dimension(0:53)             :: do_move, cube
  106.         integer, intent(in)                  :: move
  107.         integer                              :: repetition
  108.        
  109.         do_move = cube_state
  110.        
  111.        
  112.         repetition = 0
  113.         if(move > 5)   repetition = 1
  114.         if(move > 11)  repetition = 2
  115.         if(move > 17)  repetition = 3  !stupide, c'est un tour complet (juste pour test)
  116.  
  117.         do m = 0, repetition
  118.         cube = do_move         
  119.  
  120.         !tourne les cubes de la face en question
  121.         do_move = turn_face(do_move, mod(move,6))  
  122.        
  123.         !tourne ceux de la couronne
  124.         select case(mod(move,6))
  125.         case(0)
  126.             do_move(15) = cube(53)
  127.             do_move(16) = cube(50)
  128.             do_move(17) = cube(47)
  129.             do_move(18) = cube(15)
  130.             do_move(21) = cube(16)
  131.             do_move(24) = cube(17)
  132.             do_move(36) = cube(24)
  133.             do_move(37) = cube(21)
  134.             do_move(38) = cube(18)
  135.             do_move(47) = cube(36)
  136.             do_move(50) = cube(37)
  137.             do_move(53) = cube(38)
  138.         case(1)
  139.             do_move(0) = cube(18)
  140.             do_move(1) = cube(19)
  141.             do_move(2) = cube(20)
  142.             do_move(18) = cube(27)
  143.             do_move(19) = cube(28)
  144.             do_move(20) = cube(29)
  145.             do_move(27) = cube(45)
  146.             do_move(28) = cube(46)
  147.             do_move(29) = cube(47)
  148.             do_move(45) = cube(0)
  149.             do_move(46) = cube(1)
  150.             do_move(47) = cube(2)
  151.         case(2)
  152.             do_move(2) = cube(38)
  153.             do_move(5) = cube(41)
  154.             do_move(8) = cube(44)
  155.             do_move(11) = cube(2)
  156.             do_move(14) = cube(5)
  157.             do_move(17) = cube(8)
  158.             do_move(27) = cube(17)
  159.             do_move(30) = cube(14)
  160.             do_move(33) = cube(11)
  161.             do_move(38) = cube(33)
  162.             do_move(41) = cube(30)
  163.             do_move(44) = cube(27)
  164.         case(3)
  165.             do_move(9) = cube(20)
  166.             do_move(10) = cube(23)
  167.             do_move(11) = cube(26)
  168.             do_move(20) = cube(44)
  169.             do_move(23) = cube(43)
  170.             do_move(26) = cube(42)
  171.             do_move(42) = cube(45)
  172.             do_move(43) = cube(48)
  173.             do_move(44) = cube(51)
  174.             do_move(45) = cube(11)
  175.             do_move(48) = cube(10)
  176.             do_move(51) = cube(9)
  177.         case(4)
  178.             do_move(6) = cube(51)
  179.             do_move(7) = cube(52)
  180.             do_move(8) = cube(53)
  181.             do_move(24) = cube(6)
  182.             do_move(25) = cube(7)
  183.             do_move(26) = cube(8)
  184.             do_move(33) = cube(24)
  185.             do_move(34) = cube(25)
  186.             do_move(35) = cube(26)
  187.             do_move(51) = cube(33)
  188.             do_move(52) = cube(34)
  189.             do_move(53) = cube(35)
  190.         case(5)
  191.             do_move(0) = cube(9)
  192.             do_move(3) = cube(12)
  193.             do_move(6) = cube(15)
  194.             do_move(9) = cube(35)
  195.             do_move(12) = cube(32)
  196.             do_move(15) = cube(29)
  197.             do_move(29) = cube(42)
  198.             do_move(32) = cube(39)
  199.             do_move(35) = cube(36)
  200.             do_move(36) = cube(0)
  201.             do_move(39) = cube(3)
  202.             do_move(42) = cube(6)
  203.         case default
  204.             print*, "move error"
  205.         end select
  206.  
  207.         end do
  208.         return
  209.     end function do_move
  210.  
  211.     function turn_face(cube, face)
  212.        
  213.         integer, dimension(0:53), intent(in) :: cube
  214.         integer, dimension(0:53)             :: turn_face
  215.         integer, intent(in)                  :: face
  216.         integer                              :: offset
  217.        
  218.         turn_face = cube
  219.         offset = 9*face
  220.  
  221.         turn_face(mod(offset,54)) = cube(mod(offset+6,54))
  222.         turn_face(mod(offset+1,54)) = cube(mod(offset+3,54))
  223.         turn_face(mod(offset+2,54)) = cube(mod(offset+0,54))
  224.         turn_face(mod(offset+3,54)) = cube(mod(offset+7,54))
  225.         turn_face(mod(offset+5,54)) = cube(mod(offset+1,54))
  226.         turn_face(mod(offset+6,54)) = cube(mod(offset+8,54))
  227.         turn_face(mod(offset+7,54)) = cube(mod(offset+5,54))
  228.         turn_face(mod(offset+8,54)) = cube(mod(offset+2,54))
  229.  
  230.     end function turn_face
  231.  
  232.     subroutine print_cube(cube)
  233.         integer, intent(in), dimension(0:53)  :: cube
  234.         integer                               :: offset
  235.  
  236.         offset = 0
  237.         do m = 0, 5
  238.             print*, cube(0+offset:8+offset)
  239.             offset = offset+9
  240.         end do
  241.  
  242.     end subroutine print_cube
  243.  
  244.     function compare(cube1, cube2)
  245.         integer, dimension(0:53), intent(in) :: cube1, cube2
  246.         integer                              :: compare
  247.            
  248.         compare = 0
  249.         do m = 0, 53
  250.             if(cube1(m) == cube2(m)) compare = compare +1
  251.             if(cube2(m) > 5) then
  252.                 if (cube1(m) < 6) compare = compare+1  !sticker = 6=> couleur pas importante
  253.             end if
  254.         end do
  255.     end function compare
  256.  
  257.     function do_sequence(cube, move_sequence, sequence_size)
  258.         integer, intent(in)                  :: sequence_size
  259.         integer, dimension(0:53), intent(in) :: cube
  260.         integer, dimension(sequence_size), intent(in)    :: move_sequence
  261.         integer, dimension(0:53)             :: do_sequence
  262.        
  263.         do_sequence = do_move(cube, move_sequence(1)) !initialisation
  264.        
  265.         if(sequence_size > 1) then
  266.             do n = 2, sequence_size
  267.                 do_sequence = do_move(do_sequence, move_sequence(n))
  268.             end do
  269.         end if
  270.        
  271.     end function do_sequence
  272.  
  273.         function entropy(cube)
  274.         !calcule l'entropie du cube, définie comme (somme(CEP) + somme(ECeP))
  275.         !CEP = corner-edge pairs
  276.         !ECeP = edges-center pairs
  277.         !Ouaiiii c'est pas la vrai entropie puisque je la maximise ;)
  278.  
  279.         integer, dimension(0:53), intent(in)    :: cube
  280.         integer                                 :: entropy
  281.        
  282.         entropy = 0
  283.        
  284.         !Corner-edges pairs
  285.         !FLU corner
  286.         if(cube(0) == cube(1) .and. cube(15)== cube(16)) entropy = entropy+1
  287.         if(cube(0) == cube(3) .and. cube(47)== cube(50)) entropy = entropy+1
  288.         if(cube(47) == cube(46) .and. cube(15) == cube(12)) entropy = entropy+1
  289.         !FRU corner
  290.         if(cube(2) == cube(5) .and. cube(18)== cube(21)) entropy = entropy+1
  291.         if(cube(2) == cube(1) .and. cube(17)== cube(16)) entropy = entropy+1
  292.         if(cube(17) == cube(14) .and. cube(18) == cube(19)) entropy = entropy+1
  293.         !FRD corner
  294.         if(cube(8) == cube(5) .and. cube(24)== cube(21)) entropy = entropy+1
  295.         if(cube(8) == cube(7) .and. cube(38)== cube(37)) entropy = entropy+1
  296.         if(cube(38) == cube(41) .and. cube(24) == cube(25)) entropy = entropy+1
  297.         !FLD corner
  298.         if(cube(6) == cube(3) .and. cube(53)== cube(50)) entropy = entropy+1
  299.         if(cube(6) == cube(7) .and. cube(36)== cube(37)) entropy = entropy+1
  300.         if(cube(36) == cube(39) .and. cube(53) == cube(52)) entropy = entropy+1
  301.         !BLU corner
  302.         if(cube(29) == cube(32) .and. cube(45)== cube(48)) entropy = entropy+1
  303.         if(cube(29) == cube(28) .and. cube(9)== cube(10)) entropy = entropy+1
  304.         if(cube(9) == cube(12) .and. cube(45) == cube(46)) entropy = entropy+1
  305.         !BRU corner
  306.         if(cube(27) == cube(28) .and. cube(11)== cube(10)) entropy = entropy+1
  307.         if(cube(27) == cube(30) .and. cube(20)== cube(23)) entropy = entropy+1
  308.         if(cube(11) == cube(14) .and. cube(20) == cube(19)) entropy = entropy+1
  309.         !BRD corner
  310.         if(cube(33) == cube(34) .and. cube(44)== cube(43)) entropy = entropy+1
  311.         if(cube(33) == cube(30) .and. cube(26)== cube(23)) entropy = entropy+1
  312.         if(cube(44) == cube(41) .and. cube(26) == cube(25)) entropy = entropy+1
  313.         !BLD corner
  314.         if(cube(35) == cube(34) .and. cube(42)== cube(43)) entropy = entropy+1
  315.         if(cube(35) == cube(32) .and. cube(51)== cube(48)) entropy = entropy+1
  316.         if(cube(42) == cube(39) .and. cube(51) == cube(52)) entropy = entropy+1
  317.  
  318.  
  319.         !Edges-centers pairs
  320.         !F
  321.         if(cube(4) == cube(1))  entropy = entropy+1
  322.         if(cube(4) == cube(3))  entropy = entropy+1
  323.         if(cube(4) == cube(5))  entropy = entropy+1
  324.         if(cube(4) == cube(7))  entropy = entropy+1
  325.         !U
  326.         if(cube(13) == cube(10))  entropy = entropy+1
  327.         if(cube(13) == cube(12))  entropy = entropy+1
  328.         if(cube(13) == cube(14))  entropy = entropy+1
  329.         if(cube(13) == cube(16))  entropy = entropy+1
  330.         !R
  331.         if(cube(22) == cube(19))  entropy = entropy+1
  332.         if(cube(22) == cube(21))  entropy = entropy+1
  333.         if(cube(22) == cube(23))  entropy = entropy+1
  334.         if(cube(22) == cube(25))  entropy = entropy+1
  335.         !B
  336.         if(cube(31) == cube(28))  entropy = entropy+1
  337.         if(cube(31) == cube(30))  entropy = entropy+1
  338.         if(cube(31) == cube(32))  entropy = entropy+1
  339.         if(cube(31) == cube(34))  entropy = entropy+1
  340.         !D
  341.         if(cube(40) == cube(37))  entropy = entropy+1
  342.         if(cube(40) == cube(39))  entropy = entropy+1
  343.         if(cube(40) == cube(41))  entropy = entropy+1
  344.         if(cube(40) == cube(43))  entropy = entropy+1
  345.         !L
  346.         if(cube(49) == cube(46))  entropy = entropy+1
  347.         if(cube(49) == cube(48))  entropy = entropy+1
  348.         if(cube(49) == cube(50))  entropy = entropy+1
  349.         if(cube(49) == cube(52))  entropy = entropy+1
  350.  
  351.     end function entropy
  352.  
  353.     function is2x2x3(cube)
  354.         !Voir commentaires fonction entropy ...
  355.         !fonction "vérifiant" l'état du 2x2x3 situé en BD en terme de CE et ECe pairs
  356.         !vaut 16 ssi le 2x2x3 est complet
  357.  
  358.         integer, dimension(0:53), intent(in)    :: cube
  359.         integer                                 :: is2x2x3
  360.        
  361.         is2x2x3 = 0
  362.        
  363.         !Corner-edges pairs
  364.         !BRD corner
  365.         if(cube(33) == cube(34) .and. cube(44)== cube(43)) is2x2x3 = is2x2x3+1
  366.         if(cube(33) == cube(30) .and. cube(26)== cube(23)) is2x2x3 = is2x2x3+1
  367.         if(cube(44) == cube(41) .and. cube(26) == cube(25)) is2x2x3 = is2x2x3+1
  368.         !BLD corner
  369.         if(cube(35) == cube(34) .and. cube(42)== cube(43)) is2x2x3 = is2x2x3+1
  370.         if(cube(35) == cube(32) .and. cube(51)== cube(48)) is2x2x3 = is2x2x3+1
  371.         if(cube(42) == cube(39) .and. cube(51) == cube(52)) is2x2x3 = is2x2x3+1
  372.  
  373.         !Edges-centers pairs
  374.         !R
  375.         !if(cube(22) == cube(19))  is2x2x3 = is2x2x3+1
  376.         !if(cube(22) == cube(21))  is2x2x3 = is2x2x3+1
  377.         if(cube(22) == cube(23))  is2x2x3 = is2x2x3+1
  378.         if(cube(22) == cube(25))  is2x2x3 = is2x2x3+1
  379.         !B
  380.         !if(cube(31) == cube(28))  is2x2x3 = is2x2x3+1
  381.         if(cube(31) == cube(30))  is2x2x3 = is2x2x3+1
  382.         if(cube(31) == cube(32))  is2x2x3 = is2x2x3+1
  383.         if(cube(31) == cube(34))  is2x2x3 = is2x2x3+1
  384.         !D
  385.         !if(cube(40) == cube(37))  is2x2x3 = is2x2x3+1
  386.         if(cube(40) == cube(39))  is2x2x3 = is2x2x3+1
  387.         if(cube(40) == cube(41))  is2x2x3 = is2x2x3+1
  388.         if(cube(40) == cube(43))  is2x2x3 = is2x2x3+1
  389.         !L
  390.         !if(cube(49) == cube(46))  is2x2x3 = is2x2x3+1
  391.         if(cube(49) == cube(48))  is2x2x3 = is2x2x3+1
  392.         !if(cube(49) == cube(50))  is2x2x3 = is2x2x3+1
  393.         if(cube(49) == cube(52))  is2x2x3 = is2x2x3+1
  394.  
  395.     end function is2x2x3
  396.  
  397.     function is2gen(cube)
  398.         !checks if a given cube is 2-gen solvable, according to two rules :
  399.         !   correct edges orientation
  400.         !   solvable position of corners (according to Sebastian Dumitrescu's demonstration)
  401.         !   two different corners linkages must be computed !
  402.         !   returns 1 if OK, 0 otherwise
  403.  
  404.         !The given 2 faces are F & U
  405.  
  406.         use dflib
  407.  
  408.         integer, dimension(0:53), intent(in)    :: cube
  409.         integer                                 :: is2gen
  410.         integer, dimension(6)                   :: corners_pos
  411.         integer, dimension(6)                   :: corner_temp
  412.         integer                                 :: edges_ori, co
  413.  
  414.         is2gen = 0
  415.         edges_ori = 0
  416.  
  417.         !checks the orientation of edges
  418.         if(cube(1) == cube(4) .or. cube(16) == cube(13)) edges_ori = edges_ori+1
  419.         if(cube(3) == cube(4) .or. cube(50) == cube(13)) edges_ori = edges_ori+1
  420.         if(cube(5) == cube(4) .or. cube(21) == cube(13)) edges_ori = edges_ori+1
  421.         if(cube(7) == cube(4) .or. cube(37) == cube(13)) edges_ori = edges_ori+1
  422.         if(cube(28) == cube(4) .or. cube(10) == cube(13)) edges_ori = edges_ori+1
  423.         if(cube(46) == cube(4) .or. cube(13) == cube(13)) edges_ori = edges_ori+1
  424.         if(cube(19) == cube(4) .or. cube(14) == cube(13)) edges_ori = edges_ori+1
  425.        
  426.         if (edges_ori < 7) return
  427.  
  428.         !corners :
  429.         !5  6  U face
  430.         !3  4
  431.         !1  2  F face
  432.         !checks the position of the corners: 1st test (link 5-1, 6-2, 3-4 and observe the pattern)
  433.         corner_temp(1) = 2**cube(6)+ 2**cube(36)+2**cube(53) !FDL positioned corner of cube
  434.         corner_temp(2) = 2**cube(8)+2**cube(38)+2**cube(24)  !FDR
  435.         corner_temp(3) = 2**cube(0)+2**cube(15)+2**cube(47)  !FUL
  436.         corner_temp(4) = 2**cube(2)+2**cube(17)+2**cube(18)  !FUR      
  437.         corner_temp(5) = 2**cube(29)+2**cube(9)+2**cube(45)  !BUL
  438.         corner_temp(6) = 2**cube(27)+2**cube(11)+2**cube(20) !BUR
  439.        
  440.         do co = 1, 6
  441.             if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(49)) corners_pos(co)=0
  442.             if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(22)) corners_pos(co)=1
  443.             if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(49)) corners_pos(co)=2
  444.             if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(22)) corners_pos(co)=2
  445.             if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(49)) corners_pos(co)=0
  446.             if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(22)) corners_pos(co)=1
  447.         end do
  448.  
  449.         !print*, "corners positions", corners_pos
  450.        
  451.         !on doit passer en scalaire pour le select case
  452.         co = 100000*corners_pos(1)+10000*corners_pos(2)+1000*corners_pos(3)+100*corners_pos(4)+10*corners_pos(5)+corners_pos(6)
  453.         select case(co)
  454.             case(12201,1221,11022,12120,10212,102210,110220,100122,102021,101202,210021,221001,211200,210102,212010,21102,2112,22011,21210,20121,120012,112002,122100,120201,121020,201120,220110,200211,201012,202101)
  455.                 is2gen = 1
  456.             case default
  457.                 return
  458.         end select
  459.  
  460.         !2nd test (link 5-2, 6-4, 3-1 and observe the pattern)
  461.                 do co = 1, 6
  462.             if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(49)) corners_pos(co)=0
  463.             if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(22)) corners_pos(co)=1
  464.             if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(49)) corners_pos(co)=0
  465.             if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(22)) corners_pos(co)=2
  466.             if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(49)) corners_pos(co)=1
  467.             if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(22)) corners_pos(co)=2
  468.         end do
  469.  
  470.         !print*, "corners positions", corners_pos
  471.        
  472.         !on doit passer en scalaire pour le select case
  473.         co = 100000*corners_pos(1)+10000*corners_pos(2)+1000*corners_pos(3)+100*corners_pos(4)+10*corners_pos(5)+corners_pos(6)
  474.         select case(co)
  475.             case(12201,1221,11022,12120,10212,102210,110220,100122,102021,101202,210021,221001,211200,210102,212010,21102,2112,22011,21210,20121,120012,112002,122100,120201,121020,201120,220110,200211,201012,202101)
  476.                 is2gen = 1
  477.             case default
  478.                 is2gen = 0
  479.         end select
  480.  
  481.    
  482.     end function is2gen
  483.  
  484.     function fobj(cube, cube_solved, moves_sequence, sequence_size, type_func )
  485.         !La fonction objectif est définie comme nmax - x, où nmax est le nombre max de stickers OK
  486.         !en parcourant la séquence de moves et x le nombre de mouvements après lesquels on y arrive
  487.         !Exemple d'appel :
  488.         !print*, fobj(cube_scrambled, cube_solved, moves_sequence, size(moves_sequence),1)
  489.  
  490.         integer, intent(in)                             :: sequence_size, type_func
  491.         integer, dimension(0:53), intent(in)            :: cube, cube_solved
  492.         integer, dimension(sequence_size), intent(in)   :: moves_sequence
  493.         integer, dimension(2)                           :: fobj
  494.         integer     :: nmax, x, ni, i
  495.         integer, dimension(0:53)                        :: cube_temp
  496.        
  497.         cube_temp = cube
  498.        
  499.         x = 0
  500.         nmax = 0                                    !approche par paires
  501.  
  502.         !nmax = compare(cube_temp, cube_solved)  !approche par stickers
  503.         if (sequence_size > 1) then
  504.             do i = 1, sequence_size
  505.                 cube_temp = do_move(cube_temp, moves_sequence(i))
  506.  
  507.                 select case(type_func)
  508.                 case (1)
  509.                     ni = compare(cube_temp, cube_solved)            !Approche par stickers
  510.                 case (2)
  511.                     ni = entropy(cube_temp)                         !Approche par paires
  512.                 case (3)
  513.                     ni = is2x2x3(cube_temp)                         !Approche par 2x2x3 en BD
  514.                 case (4)
  515.                     ni = 10*is2gen(cube_temp) + 10*is2x2x3(cube_temp)   !Phase getting into 2-gen
  516.                 end select
  517.  
  518.                 if (ni > nmax) then
  519.                     nmax = ni
  520.                     x = i
  521.                 end if
  522.                 !print*, i ,ni, nmax, x
  523.             end do
  524.         end if
  525.  
  526.         fobj(1) = 10*nmax -x
  527.         fobj(2) = x
  528.  
  529.     end function fobj
  530.  
  531.     subroutine solve_twogen_bourrin(cube_2gen, cube_solved)
  532.     !résoud un cube dans le groupe twogen de manière optimale (HTM)
  533.     integer, dimension(0:53), intent(in)    :: cube_2gen, cube_solved
  534.     integer, dimension(0:24)                :: sol
  535.     integer, dimension(3)                   :: m
  536.     integer, dimension(2)                   :: test_resolu
  537.     integer                                 :: m1,m2,m3
  538.    
  539.     sol=sol+3   !On remplit la sol du 2-gen avec des mouvements qui n'ont rien à y faire (initialisation
  540.  
  541.  
  542.     do m1=0,5
  543.         sol(0)=(m1/2)*6 + mod(m1,2)  !0,1,6,7,12,13 OK
  544.         !print*, mod(m1,2), " sol(0) : ", sol(0)
  545.         test_resolu = fobj(cube_2gen, cube_solved, sol, 25, 1)
  546.         !print*, (test_resolu(1)+test_resolu(2))/10.
  547.         if((test_resolu(1)+test_resolu(2))/10. == 54) then
  548.             print*, "youpie 1 !"
  549.             print*, sol
  550.             GOTO 3001
  551.         end if
  552.        
  553.         do m2=0,5 !C'est pas top, ainsi il ajoute  de toute façon #m mouvements, changer les boucles...
  554.             sol(1)=(m2/2)*6 + mod(m2,2)  !0,1,6,7,12,13 OK
  555.             test_resolu = fobj(cube_2gen, cube_solved, sol, 25, 1)
  556.             if((test_resolu(1)+test_resolu(2))/10. == 54) then
  557.                 print*, "youpie 2 !"
  558.                 print*, sol
  559.                 GOTO 3001
  560.             end if
  561.         end do
  562.    
  563.     end do
  564.  
  565.    
  566. 3001    end subroutine solve_twogen_bourrin
  567.  
  568.  
  569.  
  570.     function classement_population(pop, val_fobj)
  571.     !classe une population selon les valeurs décroissantes de la fonction objectif
  572.    
  573.     integer, dimension(:, :), intent(in)    :: pop
  574.     integer, dimension(:)                   :: val_fobj   !seulement la valeur de fobj
  575.     integer                                 :: max_indice, taille, i
  576.     integer, dimension(size(pop, 1), size(pop,2))   :: classement_population
  577.  
  578.     taille = size(pop, 1)
  579.  
  580.     classement_population = pop
  581.  
  582.     do i = 1, taille
  583.         max_indice = maxloc(val_fobj,1)
  584.         classement_population(i,:) = pop(max_indice, :)
  585.         val_fobj(max_indice) = -10000
  586.     end do
  587.  
  588.     end function classement_population
  589.  
  590.     function pop_mariee (pop, phi, pc, perfo)
  591. !Choisit des couples parmi une population d'individus triés par ordre de fobj croissante:
  592. !pc = probabilité de croisement
  593. !n = taille population
  594. !phi = nbre moyen de descendants de l'Elite
  595.  
  596.     use dflib
  597.     use dfport
  598.  
  599.     integer, dimension(:,:), intent(in)             :: pop
  600.     integer, dimension(size(pop,1),2), intent(in)   :: perfo
  601.     real, intent(in)                                :: phi, pc
  602.     integer, dimension(size(pop,1),size(pop,2))     :: pop_mariee
  603.     real, dimension(0:size(pop,1)-1)                :: pi
  604.     integer                                         :: i, j, n, numero_enf, r0, p1, p2, cross
  605.     real                                            :: x0, y0
  606.    
  607.     call seed(RND$TIMESEED)
  608.     n = size(pop, 1)
  609.  
  610.     !pop2 se remplira avec les enfants
  611.     pop_mariee = pop
  612.  
  613.     !Calcul de la proba du nombre d'enfants en fct du rang
  614.  
  615.     do i = 0, n-1
  616.  
  617.         pi(i) = 1./n * (phi - i*(2.*phi-2)/(n-1.))    
  618.     end do
  619.  
  620.     !Stratégie élitiste : on conserve l'Elite -> boucle start @ numero_enf = 2 et génère (n-1)/2 couples
  621.  
  622.     numero_enf = 2
  623.     do j = 1, (n-1)/2
  624.  
  625.         !Choix du premier parent :
  626.     1   call random(x0)
  627.         r0 = int(n*0.9999*x0)            !le 0.9999 c'est pour si jamais x0==1
  628.         call random(y0)
  629.  
  630.         if( (2. - phi + 2.*y0*(phi-1.))/n > pi(r0)) then
  631.             GOTO 1
  632.         end if
  633.  
  634.         p1 = r0+1
  635.  
  636.         !Choix du 2e parent : aléatoire équiprobable
  637.  
  638.         call random(x0)
  639.         p2=int(n*0.9999*x0)+1
  640.  
  641.         !print*, p1, p2
  642.  
  643.         !Choix du gène de cross-over appartient à {1, ..., long_code-1}
  644.         !ALEATOIRE
  645.         call random(x0)
  646.         cross = (size(pop,2)-1) * 0.9999* x0 + 1   !pour le cas foireux ou x0 ==1 pile
  647.         !OU : APRES LE 1ER MAX DE FOBJ pour le premier parent
  648.         !cross = perfo(p1, 2)
  649.  
  650.         if (cross == size(pop,2)) cross = cross-1
  651.         !On attribue les gênes 1:cross du premier parent à l'individu numero_enf et cross+1:fin du 2e parent
  652.         !vice-versa à l'individu suivant (son frère)
  653.         !croisement a lieu avec une proba pc
  654.  
  655.         call random(x0)
  656.  
  657.         if(x0 < pc) then
  658.         pop_mariee(numero_enf, 1:cross) = pop(p1, 1:cross)
  659.         pop_mariee(numero_enf, cross+1:size(pop,2)) = pop(p2, cross+1:size(pop,2))
  660.         pop_mariee(numero_enf+1, 1:cross) = pop(p2, 1:cross)
  661.         pop_mariee(numero_enf+1, cross+1:size(pop,2)) = pop(p1, cross+1:size(pop,2))
  662.        
  663.         else  !pas croisement
  664.         pop_mariee(numero_enf, :)   = pop(p1, :)
  665.         pop_mariee(numero_enf+1, :) = pop(p2, :)
  666.  
  667.         end if
  668.  
  669.         numero_enf =numero_enf + 2
  670.  
  671.     end do   !boucle j sur le remplissage de pop2
  672.  
  673.     end function pop_mariee
  674.  
  675.     function pop_mutee(pop, pm)
  676.  
  677.         use dflib
  678.         use dfport
  679.  
  680.         integer, dimension(:,:), intent(in)             :: pop
  681.         real, intent(in)                                :: pm
  682.         integer, dimension(size(pop,1),size(pop,2))     :: pop_mutee
  683.         integer                                         :: i, j, n1, n2, g0
  684.         real                                            :: x0
  685.        
  686.         call seed(RND$TIMESEED)
  687.         n1 = size(pop, 1)
  688.         n2 = size(pop, 2)
  689.  
  690.         pop_mutee = pop
  691.  
  692.         do i = 2, n1  !elite ne mute pas !
  693.          do j = 1, n2
  694.             call random(x0)
  695.             if (x0 < pm) then
  696.                 call random(x0)
  697.                 g0 = 18*0.9999*x0
  698.                 pop_mutee(i,j) = g0
  699.             end if
  700.          end do
  701.         end do
  702.  
  703.  
  704.     end function pop_mutee
  705.  
  706.     function trim_sequence (moves_sequence)
  707.  
  708.         use dflib
  709.         use dfport
  710.  
  711.         integer, dimension(:), intent(in)           :: moves_sequence
  712.         integer, dimension(size(moves_sequence))    :: trim_sequence
  713.         integer :: l, ltot, shoots, move, move_suivant
  714.         real    :: r
  715.        
  716.         call seed(RND$TIMESEED)
  717.  
  718.         ltot = size(moves_sequence)
  719.         trim_sequence = moves_sequence
  720.  
  721.  
  722.         shoots = 0
  723. 21      do l = 1, ltot-1
  724.             move = trim_sequence(l)
  725.             move_suivant = trim_sequence(l+1)
  726.             if(mod(move, 6) == mod(move_suivant, 6)) then
  727.                 if(abs(move-move_suivant)== 12) then !Deux moves s'annulent
  728.                     trim_sequence(l:ltot-2) = trim_sequence(l+2:ltot)
  729.                     do shoots = 1,2
  730.                     call random(r)
  731.                     trim_sequence(ltot-2+shoots) = 18 * 0.9999* r
  732.                     end do
  733.                 end if
  734.  
  735.                 if(move < 6) then  !Face
  736.                     if(move_suivant == move) then ! face puis face = face^2
  737.                         trim_sequence(l) = move+6
  738.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  739.                         call random(r)
  740.                         trim_sequence(ltot) = 18 * 0.9999* r
  741.                     else if(move_suivant == move+6) then ! face puis face^2 = face'                
  742.                         trim_sequence(l) = move+12
  743.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  744.                         call random(r)
  745.                         trim_sequence(ltot) = 18 * 0.9999* r
  746.                     end if
  747.                 end if
  748.                
  749.                 if(move < 12) then
  750.                     if(move<6) GOTO 22
  751.                     if (move_suivant == move-6) then !Face^2 puis face = face'
  752.                         trim_sequence(l) = move+6
  753.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  754.                         call random(r)
  755.                         trim_sequence(ltot) = 18 * 0.9999* r
  756.                     else if (move_suivant == move) then !Face^2 puis face^2 = annule
  757.                         trim_sequence(l:ltot-2) = trim_sequence(l+2:ltot)
  758.                         do shoots = 1,2
  759.                         call random(r)
  760.                         trim_sequence(ltot-2+shoots) = 18 * 0.9999* r
  761.                         end do
  762.                     else if (move_suivant == move+6) then !Face^2 puis face' = face
  763.                         trim_sequence(l) = move-6
  764.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  765.                         call random(r)
  766.                         trim_sequence(ltot) = 18 * 0.9999* r
  767.                     end if
  768.                 end if
  769.                 if (move < 12)  GOTO 22
  770.                 ! Arrive ci-dessous si et seulement move et face'
  771.                 if(move_suivant == move-6)  then !face' puis face^2 = face
  772.                         trim_sequence(l) = move-12
  773.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  774.                         call random(r)
  775.                         trim_sequence(ltot) = 18 * 0.9999* r
  776.                 else if(move_suivant == move) then !face' puis face' = face^2
  777.                         trim_sequence(l) = move-6
  778.                         trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
  779.                         call random(r)
  780.                         trim_sequence(ltot) = 18 * 0.9999* r
  781.                 end if
  782. 22          GO TO 21  ! Juste pour avoir un GO TO
  783.             end if
  784.  
  785.         end do
  786.  
  787. end function trim_sequence
  788.  
  789.     function twist_sequence(moves_sequence)
  790.         !twists moves sequence clockwise around FUR corner
  791.         !vouaiii un peu stupide, autant twister le cube ...
  792.         integer, dimension(:), intent(in)       :: moves_sequence
  793.         integer, dimension(size(moves_sequence)):: twist_sequence
  794.         twist_sequence = 6*(moves_sequence/6) + mod(1+moves_sequence-6*(moves_sequence/6) ,3) + 3*((moves_sequence-6*(moves_sequence/6))/3)
  795.  
  796.     end function twist_sequence
  797.  
  798.     function twist_cube(c)
  799.         !twists whole cube around FUR corner
  800.         integer, dimension(0:53), intent(in)    :: c
  801.         integer, dimension(0:53)                :: twist_cube
  802.        
  803.         !print*, "FUR clockwise twist (= x y)"
  804.         twist_cube(0:8)   = [c(24),c(21),c(18),c(25),c(22),c(19),c(26),c(23),c(20)] !F<-R
  805.         twist_cube(9:17)  = [c(6),c(3),c(0),c(7),c(4),c(1),c(8),c(5),c(2)]          !U<-F
  806.         twist_cube(18:26) = [c(17),c(16),c(15),c(14),c(13),c(12),c(11),c(10),c(9)]  !R<-U
  807.         twist_cube(27:35) = [c(47),c(50),c(53),c(46),c(49),c(52),c(45),c(48),c(51)] !B<-L
  808.         twist_cube(36:44) = [c(33),c(30),c(27),c(34),c(31),c(28),c(35),c(32),c(29)] !D<-B
  809.         twist_cube(45:53) = c(36:44)                                                !L<-D
  810.  
  811.     end function twist_cube
  812.  
  813.     function rotate_whole_cube(cu)
  814.             !rotates whole cube clockwise around vertical axis
  815.         integer, dimension(0:53), intent(in)    :: cu
  816.         integer, dimension(0:53)                :: rotate_whole_cube
  817.        
  818.         !print*, "whole cube clockwise rotation around Oz (= y)"
  819.         rotate_whole_cube(9:17) = [cu(15),cu(12),cu(9),cu(16),cu(13),cu(10),cu(17),cu(14),cu(11)] !U
  820.         rotate_whole_cube(36:44)= [cu(38),cu(41),cu(44),cu(37),cu(40),cu(43),cu(36),cu(39),cu(42)]!D       
  821.         rotate_whole_cube(0:8)   = cu(18:26)        !F<-R
  822.         rotate_whole_cube(18:26) = cu(27:35)        !R<-B
  823.         rotate_whole_cube(27:35) = cu(45:53)        !B<-L
  824.         rotate_whole_cube(45:53) = cu(0:8)          !L<-F
  825.  
  826.     end function rotate_whole_cube
  827.  
  828.     function rotate_whole_cube_x(cu)
  829.             !rotates whole cube clockwise LR axis (= x)
  830.             !defined as xy yyy
  831.  
  832.         integer, dimension(0:53), intent(in)    :: cu
  833.         integer, dimension(0:53)                :: rotate_whole_cube_x
  834.         integer                                 :: rot
  835.                
  836.         !print*, "whole cube clockwise rotation around Oy (= x)"
  837.         rotate_whole_cube_x = cu
  838.  
  839.         rotate_whole_cube_x = twist_cube(cu)
  840.         do rot = 1,3
  841.             rotate_whole_cube_x = rotate_whole_cube(rotate_whole_cube_x)
  842.         end do
  843.  
  844.     end function rotate_whole_cube_x
  845.  
  846.    
  847.  
  848. end module external_fcts
  849.  
  850. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  851. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  852.     program RUBIK
  853.  
  854.  
  855.    
  856.     use external_fcts
  857.     use dflib
  858.     use dfport
  859.  
  860.     implicit none
  861.  
  862.     integer                                 :: i, j, k, ll, lll, two_gen, loc2x2x3, scramble_type, moves_total, twogen_bourrin
  863.     integer, dimension(:), allocatable      :: moves_sequence, scramble_sequence
  864.     integer, dimension(0:53)                :: cube_solved, cube_scrambled, cube, cube_scrambled_copy
  865.     real                                    :: r !random number
  866.     character, dimension(:),allocatable     :: scramble_notation, sequence_notation, sol1, sol2, sol3
  867.     character*2                             :: char_buffer
  868.     integer, dimension(12)                  :: best2x2x3loc, best2x2x3fobj, best2x2x3fobj_copie
  869.     integer, dimension(12, 30)              :: best2x2x3elits
  870.     integer, dimension(12,2)                :: best_temp
  871.    
  872.     !Paramètres algo génétique
  873.     integer                                     :: taille_pop, long_code, Ttot, generation, nbre_intrus, is_continuing
  874.     real                                        :: pm, pc, phi
  875.     integer, dimension(:,:), allocatable        :: pop, perfo, perfo_copie
  876.     integer, dimension(:,:,:), allocatable      :: pop_storage   !si on doit recommencer les 2x2x3, on ne perd pas tout ;-)
  877.  
  878.     !initialisation
  879.     call seed(RND$TIMESEED)   !la seed est effectivement différente à chaque appel :)
  880.     cube_solved = init(cube)
  881.     scramble_type = 2    !0 = défini manuellement, 1 = random, 2 = entré manuellement
  882.    
  883.     print*, 'Genetic rubik''s cube solver, version 2.1'
  884.     print*, ''
  885.     print*, 'Cyril Castella, May 22, 2005'
  886.     print*, 'http://www.francocube.com'
  887.     print*, ''
  888.     print*, 'Enter scramble type : 1 = random, 2 = User-defined'
  889.     read*, scramble_type
  890.     print*, ''
  891.     print*, 'Enter number of generations for simulation : '
  892.     read*, Ttot
  893.  
  894.     two_gen = 0
  895.     twogen_bourrin=0        !résoud le 2-gen par brute force au lieu algo génétique
  896.  
  897.     open(unit=10, file='results.txt', status='replace')
  898.  
  899.     !mélange
  900.  
  901.     allocate(scramble_sequence(37))
  902.     !scramble_sequence = [6,1,14,17,0,1,14,17,0,1,14,17,0,1,14,17,0,9,13,14,10,0,2,9,7,0,2,9,7,0,2,9,7,0,2,9,7,4] !FMC#81, 38 moves
  903.     !scramble_sequence = [14,16,14,6,3,7,17,16,0,7,17,16,0,7,17,16,0,8,1,5,0,7,14,5,9,12,2,10,1,12,2,10,1,12,2,10,1,11,12,15,4] !41 moves
  904.     !scramble_sequence = [0,1,6,13,0,1,6,13,0,1,6,13,6,7,12,1,6,13,0,1,6,13,0,1,6,13,0,7,12,1,6,13,12,7,0,13,12,7,0,13,12,7,0] !43 moves 2-.gen
  905.     !scramble_sequence = [10,14,9,10,17,6,15,10,12,17,6,15,10,12,17,6,15,10,12,13,3,7,15,12,2,12,3,1,10,2,12,3,1,10,2,12,3,1,10,0,1,15,12] !#82, 43 moves
  906.     scramble_sequence = [2,4,2,3,11,4,15,5,0,4,15,5,0,4,15,5,0,8,4,12,7,3,2,16,7,3,2,16,7,3,2,16,11,2,4,7,6] !FMC#83, 37 moves
  907.     !if(two_gen) scramble_sequence = mod(scramble_sequence,2) + 6*(scramble_sequence/6) !2-gen scramble
  908.  
  909.     if(scramble_type == 1) then
  910.         print*, 'Enter number of moves of the scramble'
  911.         read*, j
  912.         deallocate(scramble_sequence)
  913.         allocate(scramble_sequence(j))
  914.         do i = 1, j
  915.             call random(r)
  916.             scramble_sequence(i) = 18*0.9999*r
  917.         end do
  918.         scramble_sequence = trim_sequence(scramble_sequence)
  919.     endif
  920.  
  921.     if(scramble_type ==2) then
  922.         print*, 'Enter number of moves of the scramble'
  923.         read*, j
  924.         print*, 'Enter the scramble moves. Only one move (HTM) per line !'
  925.         deallocate(scramble_sequence)
  926.         allocate(scramble_sequence(j))
  927.         do i = 1, j
  928.             read*, char_buffer
  929.             call ToNumbers(char_buffer, scramble_sequence(i))
  930.         end do
  931.     end if
  932.  
  933.     allocate(scramble_notation(3*size(scramble_sequence)))
  934.     call ToString(scramble_sequence, scramble_notation)
  935.     print*, scramble_notation
  936.  
  937.     !Ecriture du scramble dans le fichier de résultats
  938.     write(10, *), 'Scramble : '
  939.     write(10, *), scramble_notation
  940.  
  941.     cube_scrambled = do_sequence(cube_solved, scramble_sequence, size(scramble_sequence))
  942.     cube_scrambled_copy = cube_scrambled
  943.  
  944.  
  945. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving 2x2x3 !!!!!!!!!!!!!!!!!!!!!!!!!!!
  946.  
  947.     is_continuing = 0
  948.     taille_pop = 51     !# individus
  949.     nbre_intrus = 5     !#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
  950.     long_code  = 30     !nombre de gènes  
  951.     !Ttot = 50      !nombre de générations
  952.     pm = 0.05           !proba mutation
  953.     pc = 0.8            !proba croisement
  954.     phi = 1.5           !poids élite
  955.  
  956.     allocate(pop_storage(12, taille_pop, long_code))
  957.     allocate(pop(taille_pop, long_code))
  958.     allocate(perfo(taille_pop, 2))
  959.     allocate(perfo_copie(taille_pop,2))
  960.  
  961. 101 do loc2x2x3 = 0,11
  962.     print*,''
  963.     print*, '------------2x2x3 position is : ', loc2x2x3, '       ----------'
  964.  
  965.     cube_scrambled = cube_scrambled_copy   !subtil ... je crois qu'une de mes fcts change cube_scrambled
  966.  
  967.     !Subtil ... mais visite effectivement toutes les 12 locations possibles ;-)
  968.     ! 0 = BD, 1 = xy = BL, 2 = x2 = FU, 3 = y = DL, 4 = yxy = FL, 5 = yx2 = UR
  969.     ! 6 = y2 = FD, 7 = y2xy = FR, 8 = y2x2 = UB, 9 = y3 = RD, 10 = y3xy = RB, 11 = y3x2UL
  970.  
  971.     do ll=1, loc2x2x3/3
  972.         cube_scrambled = rotate_whole_cube(cube_scrambled)
  973.         print*, "whole cube clockwise rotation around Oz (= y)"
  974.         print*
  975.     end do
  976.  
  977.     select case (modulo(loc2x2x3,3))
  978.         case(1)
  979.             cube_scrambled = twist_cube(cube_scrambled)
  980.             print*, "whole cube rotation around FUR corner (=xy)"
  981.         case(2)
  982.             do i = 1,2
  983.                 cube_scrambled = rotate_whole_cube_x(cube_scrambled)
  984.                 print*, "whole cube rotation around LR horizontal axis (=x)"
  985.             end do
  986.     end select
  987.  
  988.     k = is2gen(cube_scrambled)
  989.     !print*, "... is 2-gen ? ", k
  990.     if(k==1) then
  991.         print*, '2-gen solve'
  992.         two_gen = 1
  993.     end if
  994.    
  995.     !génération population initiale
  996.     do i = 1, taille_pop
  997.         do j = 1, long_code
  998.             call random(r)
  999.             pop(i,j) = r*0.9999*18
  1000.         end do
  1001.     end do
  1002.    
  1003.     !si on est en train de continuer car pas assez de générations, on reprend la pop d'avant pour
  1004.     !cette location
  1005.     if(is_continuing) pop = pop_storage(loc2x2x3+1, :, :)
  1006.  
  1007.     do generation = 1, Ttot
  1008.    
  1009.     if(two_gen) pop = mod(pop,2) + 6*(pop/6)        !2-gen-> séquence contient uniquement (F,U)
  1010.     !Raccourci des mouvements qui s'annulent
  1011.     do i = 1, taille_pop
  1012.         pop(i,:) = trim_sequence(pop(i,:))
  1013.     end do
  1014.    
  1015.     !Evaluation de la fonction objectif sur la population
  1016.     !pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
  1017.     do i = 1, taille_pop
  1018.     perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 3)
  1019.     !print*, pop(i,:), perfo(i, 1)
  1020.     end do
  1021.    
  1022.     !classement de la population
  1023.     perfo_copie = perfo
  1024.     pop = classement_population(pop, perfo_copie(:,1))
  1025.     perfo_copie = perfo
  1026.     perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)
  1027.  
  1028.     !Print de l'état actuel
  1029.     if (modulo(generation, 100000) == 0) then
  1030.     print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
  1031.     end if
  1032.  
  1033.     !génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
  1034.     if(nbre_intrus > 0) then
  1035.         do i = taille_pop-nbre_intrus+1, taille_pop
  1036.             do j = 1, long_code
  1037.                 call random(r)
  1038.                 pop(i,j) = r*0.9999*18
  1039.             end do
  1040.         end do
  1041.     end if
  1042.  
  1043.     !Mariages
  1044.     pop = pop_mariee (pop, phi, pc, perfo)
  1045.  
  1046.     !mutations
  1047.     pop = pop_mutee (pop, pm)
  1048.  
  1049.  
  1050.     end do  !fin boucle sur les generations
  1051.    
  1052.     pop_storage(loc2x2x3+1, :,:) = pop
  1053.  
  1054.     print*, "optimal :"
  1055.     allocate(sequence_notation(3*perfo(1,2)))
  1056.     call ToString(pop(1,1:perfo(1,2)), sequence_notation)
  1057.     print*, sequence_notation
  1058.     print*, "perfo, # moves, entropy"
  1059.     print*, perfo(1,:), (perfo(1,1)+perfo(1,2))/10
  1060.     if((perfo(1,1)+perfo(1,2))/10==16) print*, 'This 2x2x3 was solved :)'
  1061.     print*,''
  1062.    
  1063.     best2x2x3loc(loc2x2x3+1) = loc2x2x3
  1064.     best2x2x3elits(loc2x2x3+1, :) = pop(1, :)  
  1065.     best2x2x3fobj(loc2x2x3+1) = -perfo(1, 2) - (16-(perfo(1,1)+perfo(1,2))/10)*100 !pénalise une chiée si 2x2x3 pas solved
  1066.    
  1067.     deallocate(sequence_notation)
  1068.     end do !FIN BOUCLE SUR loc2x2x3
  1069.  
  1070. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Getting into 2-gen !!!!!!!!!!!!!!!!!!!!!!!!!!!
  1071.  
  1072.    
  1073.     !Classement des élites en fonction de leur fobj définie ci-dessus
  1074.  
  1075. 200 print*, ''
  1076.     best2x2x3fobj_copie = best2x2x3fobj
  1077.     best2x2x3elits = classement_population(best2x2x3elits, best2x2x3fobj_copie)
  1078.  
  1079.     best2x2x3fobj_copie = best2x2x3fobj
  1080.     best_temp(:,1) = best2x2x3loc
  1081.     best_temp(:,2) = best2x2x3fobj
  1082.     best_temp = classement_population(best_temp, best2x2x3fobj_copie)
  1083.     best2x2x3loc = best_temp(:,1)
  1084.    
  1085.     !si on n'a résolu aucun 2x2x3, ou un en plus de 20 moves on doit recommencer !!
  1086.     if(best_temp(1,2) < -20) then
  1087.         print*, 'Adding generations to solve the 2x2x3: no good solution was found.'
  1088.         is_continuing = 1
  1089.         GOTO 101
  1090.     end if
  1091.  
  1092.     !sinon, c'est cool, on continue
  1093.     deallocate(pop, perfo, perfo_copie)
  1094.    
  1095.     print*, 'Best 2x2x3 locations : '
  1096.     print*, best2x2x3loc
  1097.     print*, ''
  1098.     print*, 'End of the 2x2x3 search.'
  1099.     print*, 'Alg will now try to go into 2-gen for the best location.'
  1100.     print*, 'Location searched : ', best2x2x3loc(1)
  1101.     print*,''
  1102.  
  1103.  
  1104.     !On réinitialise cube_scrambled avec scramble, puis setup-rotations, puis best 2x2x3 solution
  1105.     loc2x2x3 = best2x2x3loc(1)
  1106.  
  1107.     cube_scrambled = cube_scrambled_copy
  1108.     do ll=1, loc2x2x3/3
  1109.         cube_scrambled = rotate_whole_cube(cube_scrambled)
  1110.         print*, "whole cube clockwise rotation around Oz (= y)"
  1111.         write(10, *), "whole cube clockwise rotation around Oz (= y)"
  1112.     end do
  1113.  
  1114.     select case (modulo(loc2x2x3,3))
  1115.         case(1)
  1116.             cube_scrambled = twist_cube(cube_scrambled)
  1117.             print*, "whole cube rotation around FUR corner (=xy)"
  1118.             write(10, *), "whole cube rotation around FUR corner (=xy)"
  1119.         case(2)
  1120.             do i = 1,2
  1121.                 cube_scrambled = rotate_whole_cube_x(cube_scrambled)
  1122.                 print*, "whole cube rotation around LR horizontal axis (=x)"
  1123.                 write(10, *), "whole cube rotation around LR horizontal axis (=x)"
  1124.             end do
  1125.     end select
  1126.  
  1127.     best_temp(1,:) = fobj(cube_scrambled, cube_solved, best2x2x3elits(1,:), long_code, 3)
  1128.     cube_scrambled = do_sequence(cube_scrambled, best2x2x3elits(1,1:best_temp(1,2)), best_temp(1,2))
  1129.     moves_total = best_temp(1,2)  !comptage pour le nombre de moves de la solution finale
  1130.  
  1131.     !Print dans le fichier de résultats
  1132.     allocate(sol1(3*moves_total))
  1133.     call ToString(best2x2x3elits(1,1:best_temp(1,2)), sol1)
  1134.  
  1135.     !!!!!!!!Algo génétique pour la phase 'getting into 2-gen'
  1136.     taille_pop = 51     !# individus
  1137.     nbre_intrus = 5     !#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
  1138.     long_code  = 30     !nombre de gènes  
  1139.     Ttot = 10000        !nombre de générations
  1140.     pm = 0.05           !proba mutation
  1141.     pc = 0.8            !proba croisement
  1142.     phi = 1.5           !poids élite
  1143.     allocate(pop(taille_pop, long_code))
  1144.     allocate(perfo(taille_pop, 2))
  1145.     allocate(perfo_copie(taille_pop,2))
  1146.  
  1147.     if(is2gen(cube_scrambled)) Ttot = 10
  1148.    
  1149.     !génération population initiale
  1150.     do i = 1, taille_pop
  1151.         do j = 1, long_code
  1152.             call random(r)
  1153.             pop(i,j) = r*0.9999*18
  1154.         end do
  1155.     end do
  1156.    
  1157.  
  1158. 201 do generation = 1, Ttot
  1159.    
  1160.     if(two_gen) pop = mod(pop,2) + 6*(pop/6)        !2-gen-> séquence contient uniquement (F,U)
  1161.     !Raccourci des mouvements qui s'annulent
  1162.     do i = 1, taille_pop
  1163.         pop(i,:) = trim_sequence(pop(i,:))
  1164.     end do
  1165.    
  1166.     !Evaluation de la fonction objectif sur la population
  1167.     !pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
  1168.     do i = 1, taille_pop
  1169.     perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 4)
  1170.     !print*, pop(i,:), perfo(i, 1)
  1171.     end do
  1172.    
  1173.     !classement de la population
  1174.     perfo_copie = perfo
  1175.     pop = classement_population(pop, perfo_copie(:,1))
  1176.     perfo_copie = perfo
  1177.     perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)
  1178.  
  1179.     !Print de l'état actuel
  1180.     if (modulo(generation, 5000) == 0) then
  1181.     print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
  1182.     end if
  1183.  
  1184.     !génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
  1185.     if(nbre_intrus > 0) then
  1186.         do i = taille_pop-nbre_intrus+1, taille_pop
  1187.             do j = 1, long_code
  1188.                 call random(r)
  1189.                 pop(i,j) = r*0.9999*18
  1190.             end do
  1191.         end do
  1192.     end if
  1193.  
  1194.     !Mariages
  1195.     pop = pop_mariee (pop, phi, pc, perfo)
  1196.  
  1197.     !mutations
  1198.     pop = pop_mutee (pop, pm)
  1199.  
  1200.  
  1201.     end do  !fin boucle sur les generations
  1202.  
  1203.    
  1204.     !Etait-on déjà dans le 2-gen groupe ?
  1205.     if(perfo(1,2)==1 .and. is2gen(cube_scrambled)==1) print*, 'Cube was already in the 2-gen group.'
  1206.     !A-t-on vraiment fini cette étape ?
  1207.  
  1208.     if( sum(perfo(1,:)) == 1700) then
  1209.        
  1210.         print*, 'Getting into 2-gen phase is done :)'
  1211.  
  1212.         allocate(sequence_notation(3*perfo(1,2)))
  1213.         call ToString(pop(1,1:perfo(1,2)), sequence_notation)
  1214.  
  1215.         !Mise à jour de l'état de cube_scrambled si on n'y était pas déjà
  1216.         if(perfo(1,2) > 1) then
  1217.             cube_scrambled = do_sequence(cube_scrambled, pop(1,1:perfo(1,2)), perfo(1,2))
  1218.             moves_total = moves_total + perfo(1,2) !comptage pour #moves solution finale
  1219.             print*, "optimal :"
  1220.             print*, sequence_notation
  1221.             print*, "# moves"
  1222.             print*, perfo(1,2)
  1223.         end if 
  1224.     else
  1225.         print*, 'Adding 10''000 gen. to solve this part'
  1226.         GOTO 201
  1227.     end if
  1228.  
  1229.     print*,''
  1230.  
  1231.     !Prépare le print final dans le fichier
  1232.     allocate(sol2(size(sequence_notation)))
  1233.     sol2 = sequence_notation
  1234.     deallocate(pop, perfo, perfo_copie, sequence_notation)
  1235.  
  1236. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving 2-gen !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1237.  
  1238. if(twogen_bourrin > 0) then
  1239.     print*, "calling bourrin"
  1240.     call solve_twogen_bourrin(cube_scrambled, cube_solved)
  1241.     GOTO 9999
  1242. end if
  1243.  
  1244. 300 taille_pop = 51     !# individus
  1245.     nbre_intrus = 5     !#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
  1246.     long_code  = 30     !nombre de gènes  
  1247.     Ttot = 200000       !nombre de générations
  1248.     pm = 0.05           !proba mutation
  1249.     pc = 0.8            !proba croisement
  1250.     phi = 1.5           !poids élite
  1251.     allocate(pop(taille_pop, long_code))
  1252.     allocate(perfo(taille_pop, 2))
  1253.     allocate(perfo_copie(taille_pop,2))
  1254.  
  1255.     two_gen = 1
  1256.  
  1257.     !génération population initiale
  1258.     do i = 1, taille_pop
  1259.         do j = 1, long_code
  1260.             call random(r)
  1261.             pop(i,j) = r*0.9999*18
  1262.         end do
  1263.     end do
  1264.    
  1265.  
  1266. 301 do generation = 1, Ttot
  1267.    
  1268.     if(two_gen) pop = mod(pop,2) + 6*(pop/6)        !2-gen-> séquence contient uniquement (F,U)
  1269.     !Raccourci des mouvements qui s'annulent
  1270.     do i = 1, taille_pop
  1271.         pop(i,:) = trim_sequence(pop(i,:))
  1272.     end do
  1273.    
  1274.     !Evaluation de la fonction objectif sur la population
  1275.     !pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
  1276.     do i = 1, taille_pop
  1277.     perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 2)
  1278.     !print*, pop(i,:), perfo(i, 1)
  1279.     end do
  1280.    
  1281.     !classement de la population
  1282.     perfo_copie = perfo
  1283.     pop = classement_population(pop, perfo_copie(:,1))
  1284.     perfo_copie = perfo
  1285.     perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)
  1286.  
  1287.     !Print de l'état actuel
  1288.     if (modulo(generation, 10000) == 0) then
  1289.     print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
  1290.     end if
  1291.  
  1292.     !génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
  1293.     if(nbre_intrus > 0) then
  1294.         do i = taille_pop-nbre_intrus+1, taille_pop
  1295.             do j = 1, long_code
  1296.                 call random(r)
  1297.                 pop(i,j) = r*0.9999*18
  1298.             end do
  1299.         end do
  1300.     end if
  1301.  
  1302.     !Mariages
  1303.     pop = pop_mariee (pop, phi, pc, perfo)
  1304.  
  1305.     !mutations
  1306.     pop = pop_mutee (pop, pm)
  1307.  
  1308.  
  1309.     end do  !fin boucle sur les generations
  1310.  
  1311.    
  1312.  
  1313.     !Si on a fini, on a fini, sinon on continue ... quel poète ce cyril !
  1314.  
  1315.     if( (perfo(1,1)+perfo(1,2))/10. == 48) then
  1316.         print*, "optimal :"
  1317.         allocate(sequence_notation(3*perfo(1,2)))
  1318.         call ToString(pop(1,1:perfo(1,2)), sequence_notation)
  1319.         print*, sequence_notation
  1320.        
  1321.         print*, "fobj / # moves / entropy"
  1322.         print*, perfo(1,:), sum(perfo(1,:))/10
  1323.         moves_total = moves_total + perfo(1,2)
  1324.        
  1325.         print*, 'Whole cube was solved in ', moves_total, ' moves :)'
  1326.     else
  1327.         print*, 'Adding 200 000 more generations to solve this stage.'
  1328.         GOTO 301
  1329.     end if
  1330.     print*,''
  1331.  
  1332.     !Prépare le print dans le fichier de résultat
  1333.     allocate(sol3(size(sequence_notation)))
  1334.     sol3 = sequence_notation
  1335.     deallocate(pop, perfo, perfo_copie, sequence_notation)
  1336.  
  1337.  
  1338.  
  1339.  
  1340.  
  1341.     !Print des résultats
  1342. 999 write(10, *) 'Solution : '
  1343.     write(10, *), sol1
  1344.     write(10, *), sol2
  1345.     write(10, *), sol3
  1346.     write(10, *), 'total moves # : ', moves_total
  1347.     close(unit=10)
  1348.  
  1349.     Deallocate(scramble_sequence, scramble_notation, sol1, sol2, sol3)
  1350.    
  1351. 9999    print*, 'End of the program.'
  1352.         read*,i
  1353.     end program RUBIK
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!