!  RUBIK.f90 
!
!  FUNCTIONS:
!	RUBIK      - Entry point of console application.
!
!

!****************************************************************************
!
!  PROGRAM: RUBIK
!
!
!****************************************************************************

module external_fcts
implicit none

integer  :: m, n, o

contains

	subroutine ToString(sequence_numbers, numbersToString)

		integer, dimension(:), intent(in)	:: sequence_numbers
		integer								:: a, taille_sequence
		character, dimension(:),intent(out)	:: numbersToString
		taille_sequence = size(sequence_numbers)
		

		do a = 1, taille_sequence
			select case(mod(sequence_numbers(a),6))
				case(0)
					numbersToString(3*a-2) = 'F'
				case(1) 
					numbersToString(3*a-2) = 'U'
				case(2) 
					numbersToString(3*a-2) = 'R'
				case(3) 
					numbersToString(3*a-2) = 'B'
				case(4) 
					numbersToString(3*a-2) = 'D'
				case(5) 
					numbersToString(3*a-2) = 'L'
			end select
			
			numbersToString(3*a-1) = ' '
			if (sequence_numbers(a)> 5) numbersToString(3*a-1) = '2'
			if (sequence_numbers(a)>11) numbersToString(3*a-1) = "'"
			numbersToString(3*a) = ' '
		end do

	end subroutine ToString

	subroutine ToNumbers(sequence_string, StringToNumbers)
		character, dimension(2), intent(in)		:: sequence_string
		integer, intent(out)					:: StringToNumbers
		integer									:: a, b, taille_string

		do a=1, 2
			select case(sequence_string(a))
				case('F')
					StringToNumbers = 0
				case('U')
					StringToNumbers = 1
				case('R')
					StringToNumbers = 2
				case('B')
					StringToNumbers = 3
				case('D')
					StringToNumbers = 4
				case('L')
					StringToNumbers = 5
				case("2")
					StringToNumbers = StringToNumbers+6
				case("'")
					StringToNumbers = StringToNumbers+12
				!case default
				!	print*, 'Incorrect caracter. Exit and run the program again'

			end select
		end do
		!print*, StringToNumbers

	end subroutine ToNumbers

	function init(cube)
		integer, dimension(0:53), intent(in)  :: cube
		integer, dimension(0:53)			  :: init
		integer								  :: indice = 0, temp = 0
		
		do m = 0, 5
			do n = 0, 8
			init(indice) = temp
			indice = indice+1
			end do
			temp = temp+1
		end do

	return
	end function init


	function do_move(cube_state, move)
		integer, dimension(0:53), intent(in) :: cube_state
		integer, dimension(0:53)			 :: do_move, cube
		integer, intent(in)					 :: move
		integer								 :: repetition
		
		do_move = cube_state
		
		
		repetition = 0
		if(move > 5)   repetition = 1
		if(move > 11)  repetition = 2
		if(move > 17)  repetition = 3  !stupide, c'est un tour complet (juste pour test)

        do m = 0, repetition
		cube = do_move  		

		!tourne les cubes de la face en question
		do_move = turn_face(do_move, mod(move,6))   
		
		!tourne ceux de la couronne
		select case(mod(move,6))
		case(0)
			do_move(15) = cube(53)
			do_move(16) = cube(50)
			do_move(17) = cube(47)
			do_move(18) = cube(15)
			do_move(21) = cube(16)
			do_move(24) = cube(17)
			do_move(36) = cube(24)
			do_move(37) = cube(21)
			do_move(38) = cube(18)
			do_move(47) = cube(36)
			do_move(50) = cube(37)
			do_move(53) = cube(38)
		case(1)
			do_move(0) = cube(18)
			do_move(1) = cube(19)
			do_move(2) = cube(20)
			do_move(18) = cube(27)
			do_move(19) = cube(28)
			do_move(20) = cube(29)
			do_move(27) = cube(45)
			do_move(28) = cube(46)
			do_move(29) = cube(47)
			do_move(45) = cube(0)
			do_move(46) = cube(1)
			do_move(47) = cube(2)
		case(2)
			do_move(2) = cube(38)
			do_move(5) = cube(41)
			do_move(8) = cube(44)
			do_move(11) = cube(2)
			do_move(14) = cube(5)
			do_move(17) = cube(8)
			do_move(27) = cube(17)
			do_move(30) = cube(14)
			do_move(33) = cube(11)
			do_move(38) = cube(33)
			do_move(41) = cube(30)
			do_move(44) = cube(27)
		case(3)
			do_move(9) = cube(20)
			do_move(10) = cube(23)
			do_move(11) = cube(26)
			do_move(20) = cube(44)
			do_move(23) = cube(43)
			do_move(26) = cube(42)
			do_move(42) = cube(45)
			do_move(43) = cube(48)
			do_move(44) = cube(51)
			do_move(45) = cube(11)
			do_move(48) = cube(10)
			do_move(51) = cube(9)
		case(4)
			do_move(6) = cube(51)
			do_move(7) = cube(52)
			do_move(8) = cube(53)
			do_move(24) = cube(6)
			do_move(25) = cube(7)
			do_move(26) = cube(8)
			do_move(33) = cube(24)
			do_move(34) = cube(25)
			do_move(35) = cube(26)
			do_move(51) = cube(33)
			do_move(52) = cube(34)
			do_move(53) = cube(35)
		case(5)
			do_move(0) = cube(9)
			do_move(3) = cube(12)
			do_move(6) = cube(15)
			do_move(9) = cube(35)
			do_move(12) = cube(32)
			do_move(15) = cube(29)
			do_move(29) = cube(42)
			do_move(32) = cube(39)
			do_move(35) = cube(36)
			do_move(36) = cube(0)
			do_move(39) = cube(3)
			do_move(42) = cube(6)
        case default
			print*, "move error"
		end select

		end do
		return
	end function do_move

	function turn_face(cube, face)
		
		integer, dimension(0:53), intent(in) :: cube
		integer, dimension(0:53)			 :: turn_face
		integer, intent(in)					 :: face
		integer								 :: offset
		
		turn_face = cube
		offset = 9*face

		turn_face(mod(offset,54)) = cube(mod(offset+6,54))
		turn_face(mod(offset+1,54)) = cube(mod(offset+3,54))
		turn_face(mod(offset+2,54)) = cube(mod(offset+0,54))
		turn_face(mod(offset+3,54)) = cube(mod(offset+7,54))
		turn_face(mod(offset+5,54)) = cube(mod(offset+1,54))
		turn_face(mod(offset+6,54)) = cube(mod(offset+8,54))
		turn_face(mod(offset+7,54)) = cube(mod(offset+5,54))
		turn_face(mod(offset+8,54)) = cube(mod(offset+2,54))

	end function turn_face 

	subroutine print_cube(cube)
		integer, intent(in), dimension(0:53)  :: cube
		integer								  :: offset

		offset = 0
		do m = 0, 5
			print*, cube(0+offset:8+offset)
			offset = offset+9
		end do

	end subroutine print_cube

	function compare(cube1, cube2)
		integer, dimension(0:53), intent(in) :: cube1, cube2
		integer								 :: compare
			
		compare = 0
		do m = 0, 53
			if(cube1(m) == cube2(m)) compare = compare +1
			if(cube2(m) > 5) then
				if (cube1(m) < 6) compare = compare+1  !sticker = 6=> couleur pas importante
			end if
		end do
	end function compare

	function do_sequence(cube, move_sequence, sequence_size)
		integer, intent(in)                  :: sequence_size
		integer, dimension(0:53), intent(in) :: cube
		integer, dimension(sequence_size), intent(in)    :: move_sequence
		integer, dimension(0:53)			 :: do_sequence
		
		do_sequence = do_move(cube, move_sequence(1)) !initialisation
		
		if(sequence_size > 1) then
			do n = 2, sequence_size
				do_sequence = do_move(do_sequence, move_sequence(n))
			end do
		end if
		
	end function do_sequence

		function entropy(cube)
		!calcule l'entropie du cube, définie comme (somme(CEP) + somme(ECeP))
		!CEP = corner-edge pairs
		!ECeP = edges-center pairs
		!Ouaiiii c'est pas la vrai entropie puisque je la maximise ;)

		integer, dimension(0:53), intent(in)	:: cube
		integer									:: entropy
		
		entropy = 0
		
		!Corner-edges pairs
		!FLU corner
		if(cube(0) == cube(1) .and. cube(15)== cube(16)) entropy = entropy+1
		if(cube(0) == cube(3) .and. cube(47)== cube(50)) entropy = entropy+1
		if(cube(47) == cube(46) .and. cube(15) == cube(12)) entropy = entropy+1
		!FRU corner
		if(cube(2) == cube(5) .and. cube(18)== cube(21)) entropy = entropy+1
		if(cube(2) == cube(1) .and. cube(17)== cube(16)) entropy = entropy+1
		if(cube(17) == cube(14) .and. cube(18) == cube(19)) entropy = entropy+1
		!FRD corner
		if(cube(8) == cube(5) .and. cube(24)== cube(21)) entropy = entropy+1
		if(cube(8) == cube(7) .and. cube(38)== cube(37)) entropy = entropy+1
		if(cube(38) == cube(41) .and. cube(24) == cube(25)) entropy = entropy+1
		!FLD corner
		if(cube(6) == cube(3) .and. cube(53)== cube(50)) entropy = entropy+1
		if(cube(6) == cube(7) .and. cube(36)== cube(37)) entropy = entropy+1
		if(cube(36) == cube(39) .and. cube(53) == cube(52)) entropy = entropy+1
		!BLU corner
		if(cube(29) == cube(32) .and. cube(45)== cube(48)) entropy = entropy+1
		if(cube(29) == cube(28) .and. cube(9)== cube(10)) entropy = entropy+1
		if(cube(9) == cube(12) .and. cube(45) == cube(46)) entropy = entropy+1
		!BRU corner
		if(cube(27) == cube(28) .and. cube(11)== cube(10)) entropy = entropy+1
		if(cube(27) == cube(30) .and. cube(20)== cube(23)) entropy = entropy+1
		if(cube(11) == cube(14) .and. cube(20) == cube(19)) entropy = entropy+1
		!BRD corner
		if(cube(33) == cube(34) .and. cube(44)== cube(43)) entropy = entropy+1
		if(cube(33) == cube(30) .and. cube(26)== cube(23)) entropy = entropy+1
		if(cube(44) == cube(41) .and. cube(26) == cube(25)) entropy = entropy+1
		!BLD corner
		if(cube(35) == cube(34) .and. cube(42)== cube(43)) entropy = entropy+1
		if(cube(35) == cube(32) .and. cube(51)== cube(48)) entropy = entropy+1
		if(cube(42) == cube(39) .and. cube(51) == cube(52)) entropy = entropy+1


		!Edges-centers pairs
		!F
		if(cube(4) == cube(1))  entropy = entropy+1
		if(cube(4) == cube(3))  entropy = entropy+1
		if(cube(4) == cube(5))  entropy = entropy+1
		if(cube(4) == cube(7))  entropy = entropy+1
		!U
		if(cube(13) == cube(10))  entropy = entropy+1
		if(cube(13) == cube(12))  entropy = entropy+1
		if(cube(13) == cube(14))  entropy = entropy+1
		if(cube(13) == cube(16))  entropy = entropy+1
		!R
		if(cube(22) == cube(19))  entropy = entropy+1
		if(cube(22) == cube(21))  entropy = entropy+1
		if(cube(22) == cube(23))  entropy = entropy+1
		if(cube(22) == cube(25))  entropy = entropy+1
		!B
		if(cube(31) == cube(28))  entropy = entropy+1
		if(cube(31) == cube(30))  entropy = entropy+1
		if(cube(31) == cube(32))  entropy = entropy+1
		if(cube(31) == cube(34))  entropy = entropy+1
		!D
		if(cube(40) == cube(37))  entropy = entropy+1
		if(cube(40) == cube(39))  entropy = entropy+1
		if(cube(40) == cube(41))  entropy = entropy+1
		if(cube(40) == cube(43))  entropy = entropy+1
		!L
		if(cube(49) == cube(46))  entropy = entropy+1
		if(cube(49) == cube(48))  entropy = entropy+1
		if(cube(49) == cube(50))  entropy = entropy+1
		if(cube(49) == cube(52))  entropy = entropy+1

	end function entropy

	function is2x2x3(cube)
		!Voir commentaires fonction entropy ...
		!fonction "vérifiant" l'état du 2x2x3 situé en BD en terme de CE et ECe pairs
		!vaut 16 ssi le 2x2x3 est complet

		integer, dimension(0:53), intent(in)	:: cube
		integer									:: is2x2x3
		
		is2x2x3 = 0
		
		!Corner-edges pairs
		!BRD corner
		if(cube(33) == cube(34) .and. cube(44)== cube(43)) is2x2x3 = is2x2x3+1
		if(cube(33) == cube(30) .and. cube(26)== cube(23)) is2x2x3 = is2x2x3+1
		if(cube(44) == cube(41) .and. cube(26) == cube(25)) is2x2x3 = is2x2x3+1
		!BLD corner
		if(cube(35) == cube(34) .and. cube(42)== cube(43)) is2x2x3 = is2x2x3+1
		if(cube(35) == cube(32) .and. cube(51)== cube(48)) is2x2x3 = is2x2x3+1
		if(cube(42) == cube(39) .and. cube(51) == cube(52)) is2x2x3 = is2x2x3+1

		!Edges-centers pairs
		!R
		!if(cube(22) == cube(19))  is2x2x3 = is2x2x3+1
		!if(cube(22) == cube(21))  is2x2x3 = is2x2x3+1
		if(cube(22) == cube(23))  is2x2x3 = is2x2x3+1
		if(cube(22) == cube(25))  is2x2x3 = is2x2x3+1
		!B
		!if(cube(31) == cube(28))  is2x2x3 = is2x2x3+1
		if(cube(31) == cube(30))  is2x2x3 = is2x2x3+1
		if(cube(31) == cube(32))  is2x2x3 = is2x2x3+1
		if(cube(31) == cube(34))  is2x2x3 = is2x2x3+1
		!D
		!if(cube(40) == cube(37))  is2x2x3 = is2x2x3+1
		if(cube(40) == cube(39))  is2x2x3 = is2x2x3+1
		if(cube(40) == cube(41))  is2x2x3 = is2x2x3+1
		if(cube(40) == cube(43))  is2x2x3 = is2x2x3+1
		!L
		!if(cube(49) == cube(46))  is2x2x3 = is2x2x3+1
		if(cube(49) == cube(48))  is2x2x3 = is2x2x3+1
		!if(cube(49) == cube(50))  is2x2x3 = is2x2x3+1
		if(cube(49) == cube(52))  is2x2x3 = is2x2x3+1

	end function is2x2x3

	function is2gen(cube)
		!checks if a given cube is 2-gen solvable, according to two rules :
		!	correct edges orientation
		!	solvable position of corners (according to Sebastian Dumitrescu's demonstration)
		!	two different corners linkages must be computed !
		!	returns 1 if OK, 0 otherwise

		!The given 2 faces are F & U

		use dflib

		integer, dimension(0:53), intent(in)	:: cube
		integer									:: is2gen
		integer, dimension(6)					:: corners_pos
		integer, dimension(6)					:: corner_temp
		integer									:: edges_ori, co

		is2gen = 0
		edges_ori = 0

		!checks the orientation of edges
		if(cube(1) == cube(4) .or. cube(16) == cube(13)) edges_ori = edges_ori+1
		if(cube(3) == cube(4) .or. cube(50) == cube(13)) edges_ori = edges_ori+1
		if(cube(5) == cube(4) .or. cube(21) == cube(13)) edges_ori = edges_ori+1
		if(cube(7) == cube(4) .or. cube(37) == cube(13)) edges_ori = edges_ori+1
		if(cube(28) == cube(4) .or. cube(10) == cube(13)) edges_ori = edges_ori+1
		if(cube(46) == cube(4) .or. cube(13) == cube(13)) edges_ori = edges_ori+1
		if(cube(19) == cube(4) .or. cube(14) == cube(13)) edges_ori = edges_ori+1
		
		if (edges_ori < 7) return

		!corners :
		!5  6  U face
		!3  4
		!1  2  F face
		!checks the position of the corners: 1st test (link 5-1, 6-2, 3-4 and observe the pattern)
		corner_temp(1) = 2**cube(6)+ 2**cube(36)+2**cube(53) !FDL positioned corner of cube
		corner_temp(2) = 2**cube(8)+2**cube(38)+2**cube(24)	 !FDR
		corner_temp(3) = 2**cube(0)+2**cube(15)+2**cube(47)	 !FUL
		corner_temp(4) = 2**cube(2)+2**cube(17)+2**cube(18)	 !FUR		
		corner_temp(5) = 2**cube(29)+2**cube(9)+2**cube(45)	 !BUL
		corner_temp(6) = 2**cube(27)+2**cube(11)+2**cube(20) !BUR
		
		do co = 1, 6
			if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(49)) corners_pos(co)=0
			if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(22)) corners_pos(co)=1
			if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(49)) corners_pos(co)=2
			if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(22)) corners_pos(co)=2
			if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(49)) corners_pos(co)=0
			if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(22)) corners_pos(co)=1
		end do

		!print*, "corners positions", corners_pos
		
		!on doit passer en scalaire pour le select case
		co = 100000*corners_pos(1)+10000*corners_pos(2)+1000*corners_pos(3)+100*corners_pos(4)+10*corners_pos(5)+corners_pos(6)	
		select case(co)
			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) 
				is2gen = 1
			case default
				return
		end select

		!2nd test (link 5-2, 6-4, 3-1 and observe the pattern)
				do co = 1, 6
			if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(49)) corners_pos(co)=0
			if (corner_temp(co) == 2**cube(4)+2**cube(40)+2**cube(22)) corners_pos(co)=1
			if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(49)) corners_pos(co)=0
			if (corner_temp(co) == 2**cube(4)+2**cube(13)+2**cube(22)) corners_pos(co)=2
			if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(49)) corners_pos(co)=1
			if (corner_temp(co) == 2**cube(31)+2**cube(13)+2**cube(22)) corners_pos(co)=2
		end do

		!print*, "corners positions", corners_pos
		
		!on doit passer en scalaire pour le select case
		co = 100000*corners_pos(1)+10000*corners_pos(2)+1000*corners_pos(3)+100*corners_pos(4)+10*corners_pos(5)+corners_pos(6)	
		select case(co)
			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) 
				is2gen = 1
			case default
				is2gen = 0
		end select

	
	end function is2gen

	function fobj(cube, cube_solved, moves_sequence, sequence_size, type_func )
		!La fonction objectif est définie comme nmax - x, où nmax est le nombre max de stickers OK
		!en parcourant la séquence de moves et x le nombre de mouvements après lesquels on y arrive
		!Exemple d'appel :
		!print*, fobj(cube_scrambled, cube_solved, moves_sequence, size(moves_sequence),1) 

		integer, intent(in)								:: sequence_size, type_func
		integer, dimension(0:53), intent(in)			:: cube, cube_solved
		integer, dimension(sequence_size), intent(in)	:: moves_sequence
		integer, dimension(2)							:: fobj
		integer		:: nmax, x, ni, i
		integer, dimension(0:53)						:: cube_temp
		
		cube_temp = cube
		
		x = 0
		nmax = 0									!approche par paires

		!nmax = compare(cube_temp, cube_solved)  !approche par stickers
		if (sequence_size > 1) then
			do i = 1, sequence_size
				cube_temp = do_move(cube_temp, moves_sequence(i))

				select case(type_func)
				case (1)
					ni = compare(cube_temp, cube_solved)			!Approche par stickers
				case (2)
					ni = entropy(cube_temp)							!Approche par paires
				case (3)
					ni = is2x2x3(cube_temp)							!Approche par 2x2x3 en BD
				case (4)
					ni = 10*is2gen(cube_temp) + 10*is2x2x3(cube_temp)	!Phase getting into 2-gen
				end select

				if (ni > nmax) then
					nmax = ni
					x = i
				end if
				!print*, i ,ni, nmax, x
			end do
		end if

		fobj(1) = 10*nmax -x
		fobj(2) = x

	end function fobj

	subroutine solve_twogen_bourrin(cube_2gen, cube_solved)
	!résoud un cube dans le groupe twogen de manière optimale (HTM)
	integer, dimension(0:53), intent(in)	:: cube_2gen, cube_solved
	integer, dimension(0:24)				:: sol
	integer, dimension(3)					:: m
	integer, dimension(2)					:: test_resolu
	integer									:: m1,m2,m3
	
	sol=sol+3   !On remplit la sol du 2-gen avec des mouvements qui n'ont rien à y faire (initialisation


	do m1=0,5
		sol(0)=(m1/2)*6 + mod(m1,2)  !0,1,6,7,12,13 OK
		!print*, mod(m1,2), " sol(0) : ", sol(0)
		test_resolu = fobj(cube_2gen, cube_solved, sol, 25, 1)
		!print*, (test_resolu(1)+test_resolu(2))/10.
		if((test_resolu(1)+test_resolu(2))/10. == 54) then
			print*, "youpie 1 !"
			print*, sol
			GOTO 3001
		end if
		
		do m2=0,5 !C'est pas top, ainsi il ajoute  de toute façon #m mouvements, changer les boucles...
			sol(1)=(m2/2)*6 + mod(m2,2)  !0,1,6,7,12,13 OK
			test_resolu = fobj(cube_2gen, cube_solved, sol, 25, 1)
			if((test_resolu(1)+test_resolu(2))/10. == 54) then
				print*, "youpie 2 !"
				print*, sol
				GOTO 3001
			end if
		end do
	
	end do

	
3001	end subroutine solve_twogen_bourrin



	function classement_population(pop, val_fobj)
	!classe une population selon les valeurs décroissantes de la fonction objectif
	
	integer, dimension(:, :), intent(in)	:: pop
	integer, dimension(:)					:: val_fobj   !seulement la valeur de fobj
	integer									:: max_indice, taille, i
	integer, dimension(size(pop, 1), size(pop,2))	:: classement_population

	taille = size(pop, 1)

	classement_population = pop

	do i = 1, taille
		max_indice = maxloc(val_fobj,1)
		classement_population(i,:) = pop(max_indice, :)
		val_fobj(max_indice) = -10000
	end do

	end function classement_population

	function pop_mariee (pop, phi, pc, perfo)
!Choisit des couples parmi une population d'individus triés par ordre de fobj croissante:
!pc = probabilité de croisement
!n = taille population
!phi = nbre moyen de descendants de l'Elite

	use dflib
	use dfport

	integer, dimension(:,:), intent(in)				:: pop
	integer, dimension(size(pop,1),2), intent(in)	:: perfo
	real, intent(in)								:: phi, pc
	integer, dimension(size(pop,1),size(pop,2))		:: pop_mariee
	real, dimension(0:size(pop,1)-1)				:: pi
	integer											:: i, j, n, numero_enf, r0, p1, p2, cross
	real											:: x0, y0
	
	call seed(RND$TIMESEED)
	n = size(pop, 1)

	!pop2 se remplira avec les enfants
	pop_mariee = pop

	!Calcul de la proba du nombre d'enfants en fct du rang

	do i = 0, n-1

		pi(i) = 1./n * (phi - i*(2.*phi-2)/(n-1.))    
	end do

	!Stratégie élitiste : on conserve l'Elite -> boucle start @ numero_enf = 2 et génère (n-1)/2 couples

	numero_enf = 2
	do j = 1, (n-1)/2

		!Choix du premier parent :
	1   call random(x0)
		r0 = int(n*0.9999*x0)			 !le 0.9999 c'est pour si jamais x0==1
		call random(y0)

		if( (2. - phi + 2.*y0*(phi-1.))/n > pi(r0)) then
			GOTO 1
		end if

		p1 = r0+1

		!Choix du 2e parent : aléatoire équiprobable

		call random(x0)
		p2=int(n*0.9999*x0)+1

		!print*, p1, p2

		!Choix du gène de cross-over appartient à {1, ..., long_code-1}
		!ALEATOIRE
		call random(x0)
		cross = (size(pop,2)-1) * 0.9999* x0 + 1   !pour le cas foireux ou x0 ==1 pile
		!OU : APRES LE 1ER MAX DE FOBJ pour le premier parent
		!cross = perfo(p1, 2)

		if (cross == size(pop,2)) cross = cross-1
		!On attribue les gênes 1:cross du premier parent à l'individu numero_enf et cross+1:fin du 2e parent
		!vice-versa à l'individu suivant (son frère)
		!croisement a lieu avec une proba pc

		call random(x0)

		if(x0 < pc) then
		pop_mariee(numero_enf, 1:cross) = pop(p1, 1:cross)
		pop_mariee(numero_enf, cross+1:size(pop,2)) = pop(p2, cross+1:size(pop,2))
		pop_mariee(numero_enf+1, 1:cross) = pop(p2, 1:cross)
		pop_mariee(numero_enf+1, cross+1:size(pop,2)) = pop(p1, cross+1:size(pop,2))
		
		else  !pas croisement
		pop_mariee(numero_enf, :)   = pop(p1, :)
		pop_mariee(numero_enf+1, :) = pop(p2, :)

		end if

		numero_enf =numero_enf + 2

	end do   !boucle j sur le remplissage de pop2

	end function pop_mariee

	function pop_mutee(pop, pm)

		use dflib
		use dfport

		integer, dimension(:,:), intent(in)				:: pop
		real, intent(in)								:: pm
		integer, dimension(size(pop,1),size(pop,2))		:: pop_mutee
		integer											:: i, j, n1, n2, g0
		real											:: x0
		
		call seed(RND$TIMESEED)
		n1 = size(pop, 1)
		n2 = size(pop, 2)

		pop_mutee = pop

		do i = 2, n1  !elite ne mute pas !
		 do j = 1, n2
			call random(x0)
			if (x0 < pm) then
				call random(x0)
				g0 = 18*0.9999*x0 
				pop_mutee(i,j) = g0
			end if
		 end do
		end do


	end function pop_mutee

	function trim_sequence (moves_sequence)

		use dflib
		use dfport

		integer, dimension(:), intent(in)			:: moves_sequence
		integer, dimension(size(moves_sequence))	:: trim_sequence
		integer	:: l, ltot, shoots, move, move_suivant
		real	:: r
		
		call seed(RND$TIMESEED)

		ltot = size(moves_sequence)
		trim_sequence = moves_sequence


		shoots = 0
21		do l = 1, ltot-1
			move = trim_sequence(l)
			move_suivant = trim_sequence(l+1)
			if(mod(move, 6) == mod(move_suivant, 6)) then
				if(abs(move-move_suivant)== 12) then !Deux moves s'annulent
					trim_sequence(l:ltot-2) = trim_sequence(l+2:ltot)
					do shoots = 1,2
					call random(r)
					trim_sequence(ltot-2+shoots) = 18 * 0.9999* r
					end do
				end if

				if(move < 6) then  !Face
					if(move_suivant == move) then ! face puis face = face^2
						trim_sequence(l) = move+6
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
					else if(move_suivant == move+6) then ! face puis face^2 = face'					
						trim_sequence(l) = move+12
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
					end if
				end if
				
				if(move < 12) then
					if(move<6) GOTO 22
					if (move_suivant == move-6) then !Face^2 puis face = face'
						trim_sequence(l) = move+6
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
					else if (move_suivant == move) then !Face^2 puis face^2 = annule
						trim_sequence(l:ltot-2) = trim_sequence(l+2:ltot)
						do shoots = 1,2
						call random(r)
						trim_sequence(ltot-2+shoots) = 18 * 0.9999* r
						end do
					else if (move_suivant == move+6) then !Face^2 puis face' = face
						trim_sequence(l) = move-6
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
					end if
				end if
				if (move < 12)  GOTO 22
				! Arrive ci-dessous si et seulement move et face'
				if(move_suivant == move-6)  then !face' puis face^2 = face
						trim_sequence(l) = move-12
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
				else if(move_suivant == move) then !face' puis face' = face^2
						trim_sequence(l) = move-6
						trim_sequence(l+1:ltot-1) = trim_sequence(l+2:ltot)
						call random(r)
						trim_sequence(ltot) = 18 * 0.9999* r
				end if
22			GO TO 21  ! Juste pour avoir un GO TO
			end if

		end do

end function trim_sequence

	function twist_sequence(moves_sequence)
		!twists moves sequence clockwise around FUR corner
		!vouaiii un peu stupide, autant twister le cube ...
		integer, dimension(:), intent(in)		:: moves_sequence
		integer, dimension(size(moves_sequence)):: twist_sequence
		twist_sequence = 6*(moves_sequence/6) + mod(1+moves_sequence-6*(moves_sequence/6) ,3) + 3*((moves_sequence-6*(moves_sequence/6))/3)

	end function twist_sequence

	function twist_cube(c)
		!twists whole cube around FUR corner
		integer, dimension(0:53), intent(in)	:: c
		integer, dimension(0:53)				:: twist_cube
		
		!print*, "FUR clockwise twist (= x y)"
		twist_cube(0:8)   = [c(24),c(21),c(18),c(25),c(22),c(19),c(26),c(23),c(20)]	!F<-R
		twist_cube(9:17)  = [c(6),c(3),c(0),c(7),c(4),c(1),c(8),c(5),c(2)]			!U<-F
		twist_cube(18:26) = [c(17),c(16),c(15),c(14),c(13),c(12),c(11),c(10),c(9)]	!R<-U
		twist_cube(27:35) = [c(47),c(50),c(53),c(46),c(49),c(52),c(45),c(48),c(51)]	!B<-L
		twist_cube(36:44) = [c(33),c(30),c(27),c(34),c(31),c(28),c(35),c(32),c(29)]	!D<-B
		twist_cube(45:53) = c(36:44)												!L<-D

	end function twist_cube

	function rotate_whole_cube(cu)
			!rotates whole cube clockwise around vertical axis
		integer, dimension(0:53), intent(in)	:: cu
		integer, dimension(0:53)				:: rotate_whole_cube
		
		!print*, "whole cube clockwise rotation around Oz (= y)"
		rotate_whole_cube(9:17) = [cu(15),cu(12),cu(9),cu(16),cu(13),cu(10),cu(17),cu(14),cu(11)] !U
		rotate_whole_cube(36:44)= [cu(38),cu(41),cu(44),cu(37),cu(40),cu(43),cu(36),cu(39),cu(42)]!D		
		rotate_whole_cube(0:8)   = cu(18:26)		!F<-R
		rotate_whole_cube(18:26) = cu(27:35)		!R<-B
		rotate_whole_cube(27:35) = cu(45:53)		!B<-L
		rotate_whole_cube(45:53) = cu(0:8)			!L<-F

	end function rotate_whole_cube

	function rotate_whole_cube_x(cu)
			!rotates whole cube clockwise LR axis (= x)
			!defined as xy yyy

		integer, dimension(0:53), intent(in)	:: cu
		integer, dimension(0:53)				:: rotate_whole_cube_x
		integer									:: rot
				
		!print*, "whole cube clockwise rotation around Oy (= x)"
		rotate_whole_cube_x = cu

		rotate_whole_cube_x = twist_cube(cu)
		do rot = 1,3
			rotate_whole_cube_x = rotate_whole_cube(rotate_whole_cube_x)
		end do

	end function rotate_whole_cube_x

	

end module external_fcts

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	program RUBIK


	
	use external_fcts
	use dflib
	use dfport

	implicit none

	integer									:: i, j, k, ll, lll, two_gen, loc2x2x3, scramble_type, moves_total, twogen_bourrin
	integer, dimension(:), allocatable		:: moves_sequence, scramble_sequence
	integer, dimension(0:53)				:: cube_solved, cube_scrambled, cube, cube_scrambled_copy
	real									:: r !random number
	character, dimension(:),allocatable		:: scramble_notation, sequence_notation, sol1, sol2, sol3
	character*2								:: char_buffer
	integer, dimension(12)					:: best2x2x3loc, best2x2x3fobj, best2x2x3fobj_copie
	integer, dimension(12, 30)				:: best2x2x3elits
	integer, dimension(12,2)				:: best_temp
	
	!Paramètres algo génétique
	integer										:: taille_pop, long_code, Ttot, generation, nbre_intrus, is_continuing
	real										:: pm, pc, phi
	integer, dimension(:,:), allocatable		:: pop, perfo, perfo_copie
	integer, dimension(:,:,:), allocatable		:: pop_storage   !si on doit recommencer les 2x2x3, on ne perd pas tout ;-)

	!initialisation
	call seed(RND$TIMESEED)   !la seed est effectivement différente à chaque appel :)
	cube_solved = init(cube)
	scramble_type = 2    !0 = défini manuellement, 1 = random, 2 = entré manuellement
	
	print*, 'Genetic rubik''s cube solver, version 2.1'
	print*, ''
	print*, 'Cyril Castella, May 22, 2005'
	print*, 'http://www.francocube.com'
	print*, ''
	print*, 'Enter scramble type : 1 = random, 2 = User-defined'
	read*, scramble_type
	print*, ''
	print*, 'Enter number of generations for simulation : '
	read*, Ttot

	two_gen = 0
	twogen_bourrin=0		!résoud le 2-gen par brute force au lieu algo génétique

	open(unit=10, file='results.txt', status='replace')

	!mélange

	allocate(scramble_sequence(37))
	!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
	!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
	!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
	!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
	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
	!if(two_gen) scramble_sequence = mod(scramble_sequence,2) + 6*(scramble_sequence/6) !2-gen scramble

	if(scramble_type == 1) then
		print*, 'Enter number of moves of the scramble'
		read*, j
		deallocate(scramble_sequence)
		allocate(scramble_sequence(j))
		do i = 1, j
			call random(r)
			scramble_sequence(i) = 18*0.9999*r
		end do
		scramble_sequence = trim_sequence(scramble_sequence)
	endif

	if(scramble_type ==2) then
		print*, 'Enter number of moves of the scramble'
		read*, j
		print*, 'Enter the scramble moves. Only one move (HTM) per line !'
		deallocate(scramble_sequence)
		allocate(scramble_sequence(j))
		do i = 1, j
			read*, char_buffer
			call ToNumbers(char_buffer, scramble_sequence(i))
		end do
	end if

	allocate(scramble_notation(3*size(scramble_sequence)))
	call ToString(scramble_sequence, scramble_notation)
	print*, scramble_notation

	!Ecriture du scramble dans le fichier de résultats
	write(10, *), 'Scramble : '
	write(10, *), scramble_notation

	cube_scrambled = do_sequence(cube_solved, scramble_sequence, size(scramble_sequence))
	cube_scrambled_copy = cube_scrambled


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving 2x2x3 !!!!!!!!!!!!!!!!!!!!!!!!!!!

	is_continuing = 0
	taille_pop = 51	    !# individus
	nbre_intrus = 5		!#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
	long_code  = 30		!nombre de gènes	
	!Ttot = 50		!nombre de générations
	pm = 0.05			!proba mutation
	pc = 0.8			!proba croisement
	phi = 1.5			!poids élite

	allocate(pop_storage(12, taille_pop, long_code))
	allocate(pop(taille_pop, long_code))
	allocate(perfo(taille_pop, 2))
	allocate(perfo_copie(taille_pop,2))

101	do loc2x2x3 = 0,11
	print*,''
	print*, '------------2x2x3 position is : ', loc2x2x3, '       ----------'

	cube_scrambled = cube_scrambled_copy   !subtil ... je crois qu'une de mes fcts change cube_scrambled

	!Subtil ... mais visite effectivement toutes les 12 locations possibles ;-)
	! 0 = BD, 1 = xy = BL, 2 = x2 = FU, 3 = y = DL, 4 = yxy = FL, 5 = yx2 = UR
	! 6 = y2 = FD, 7 = y2xy = FR, 8 = y2x2 = UB, 9 = y3 = RD, 10 = y3xy = RB, 11 = y3x2UL

	do ll=1, loc2x2x3/3
		cube_scrambled = rotate_whole_cube(cube_scrambled)
		print*, "whole cube clockwise rotation around Oz (= y)"
		print*
	end do

	select case (modulo(loc2x2x3,3))
		case(1)
			cube_scrambled = twist_cube(cube_scrambled)
			print*, "whole cube rotation around FUR corner (=xy)"
		case(2)
			do i = 1,2
				cube_scrambled = rotate_whole_cube_x(cube_scrambled)
				print*, "whole cube rotation around LR horizontal axis (=x)"
			end do
	end select

	k = is2gen(cube_scrambled)
	!print*, "... is 2-gen ? ", k
	if(k==1) then 
		print*, '2-gen solve'
		two_gen = 1
	end if
	
	!génération population initiale
	do i = 1, taille_pop
		do j = 1, long_code
			call random(r)
			pop(i,j) = r*0.9999*18
		end do
	end do
	
	!si on est en train de continuer car pas assez de générations, on reprend la pop d'avant pour
	!cette location
	if(is_continuing) pop = pop_storage(loc2x2x3+1, :, :)

	do generation = 1, Ttot
	
	if(two_gen) pop = mod(pop,2) + 6*(pop/6)		!2-gen-> séquence contient uniquement (F,U)
	!Raccourci des mouvements qui s'annulent
	do i = 1, taille_pop
		pop(i,:) = trim_sequence(pop(i,:))
	end do
	
	!Evaluation de la fonction objectif sur la population
	!pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
	do i = 1, taille_pop
	perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 3) 
	!print*, pop(i,:), perfo(i, 1)
	end do
	
	!classement de la population
	perfo_copie = perfo
	pop = classement_population(pop, perfo_copie(:,1))
	perfo_copie = perfo
	perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)

	!Print de l'état actuel
	if (modulo(generation, 100000) == 0) then
	print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
	end if

	!génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
	if(nbre_intrus > 0) then 
		do i = taille_pop-nbre_intrus+1, taille_pop
			do j = 1, long_code
				call random(r)
				pop(i,j) = r*0.9999*18
			end do
		end do
	end if

	!Mariages
	pop = pop_mariee (pop, phi, pc, perfo)

	!mutations
	pop = pop_mutee (pop, pm)


	end do  !fin boucle sur les generations
	
	pop_storage(loc2x2x3+1, :,:) = pop

	print*, "optimal :"
	allocate(sequence_notation(3*perfo(1,2)))
	call ToString(pop(1,1:perfo(1,2)), sequence_notation)
	print*, sequence_notation
	print*, "perfo, # moves, entropy"
	print*, perfo(1,:), (perfo(1,1)+perfo(1,2))/10
	if((perfo(1,1)+perfo(1,2))/10==16) print*, 'This 2x2x3 was solved :)'
	print*,''
	
	best2x2x3loc(loc2x2x3+1) = loc2x2x3
	best2x2x3elits(loc2x2x3+1, :) = pop(1, :)	
	best2x2x3fobj(loc2x2x3+1) = -perfo(1, 2) - (16-(perfo(1,1)+perfo(1,2))/10)*100 !pénalise une chiée si 2x2x3 pas solved
	
	deallocate(sequence_notation)
	end do !FIN BOUCLE SUR loc2x2x3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Getting into 2-gen !!!!!!!!!!!!!!!!!!!!!!!!!!!

	
	!Classement des élites en fonction de leur fobj définie ci-dessus

200	print*, ''
	best2x2x3fobj_copie = best2x2x3fobj
	best2x2x3elits = classement_population(best2x2x3elits, best2x2x3fobj_copie)

	best2x2x3fobj_copie = best2x2x3fobj
	best_temp(:,1) = best2x2x3loc
	best_temp(:,2) = best2x2x3fobj
	best_temp = classement_population(best_temp, best2x2x3fobj_copie)
	best2x2x3loc = best_temp(:,1)
	
	!si on n'a résolu aucun 2x2x3, ou un en plus de 20 moves on doit recommencer !!
	if(best_temp(1,2) < -20) then
		print*, 'Adding generations to solve the 2x2x3: no good solution was found.'
		is_continuing = 1
		GOTO 101
	end if

	!sinon, c'est cool, on continue
	deallocate(pop, perfo, perfo_copie)
	
	print*, 'Best 2x2x3 locations : '
	print*, best2x2x3loc
	print*, ''
	print*, 'End of the 2x2x3 search.'
	print*, 'Alg will now try to go into 2-gen for the best location.'
	print*, 'Location searched : ', best2x2x3loc(1)
	print*,''


	!On réinitialise cube_scrambled avec scramble, puis setup-rotations, puis best 2x2x3 solution
	loc2x2x3 = best2x2x3loc(1)

	cube_scrambled = cube_scrambled_copy
	do ll=1, loc2x2x3/3
		cube_scrambled = rotate_whole_cube(cube_scrambled)
		print*, "whole cube clockwise rotation around Oz (= y)"
		write(10, *), "whole cube clockwise rotation around Oz (= y)"
	end do

	select case (modulo(loc2x2x3,3))
		case(1)
			cube_scrambled = twist_cube(cube_scrambled)
			print*, "whole cube rotation around FUR corner (=xy)"
			write(10, *), "whole cube rotation around FUR corner (=xy)"
		case(2)
			do i = 1,2
				cube_scrambled = rotate_whole_cube_x(cube_scrambled)
				print*, "whole cube rotation around LR horizontal axis (=x)"
				write(10, *), "whole cube rotation around LR horizontal axis (=x)"
			end do
	end select

	best_temp(1,:) = fobj(cube_scrambled, cube_solved, best2x2x3elits(1,:), long_code, 3)
	cube_scrambled = do_sequence(cube_scrambled, best2x2x3elits(1,1:best_temp(1,2)), best_temp(1,2))
	moves_total = best_temp(1,2)  !comptage pour le nombre de moves de la solution finale

	!Print dans le fichier de résultats
	allocate(sol1(3*moves_total))
	call ToString(best2x2x3elits(1,1:best_temp(1,2)), sol1)

	!!!!!!!!Algo génétique pour la phase 'getting into 2-gen'
	taille_pop = 51	    !# individus
	nbre_intrus = 5		!#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
	long_code  = 30		!nombre de gènes	
	Ttot = 10000		!nombre de générations
	pm = 0.05			!proba mutation
	pc = 0.8			!proba croisement
	phi = 1.5			!poids élite
	allocate(pop(taille_pop, long_code))
	allocate(perfo(taille_pop, 2))
	allocate(perfo_copie(taille_pop,2))

	if(is2gen(cube_scrambled)) Ttot = 10
	
	!génération population initiale
	do i = 1, taille_pop
		do j = 1, long_code
			call random(r)
			pop(i,j) = r*0.9999*18
		end do
	end do
	

201	do generation = 1, Ttot
	
	if(two_gen) pop = mod(pop,2) + 6*(pop/6)		!2-gen-> séquence contient uniquement (F,U)
	!Raccourci des mouvements qui s'annulent
	do i = 1, taille_pop
		pop(i,:) = trim_sequence(pop(i,:))
	end do
	
	!Evaluation de la fonction objectif sur la population
	!pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
	do i = 1, taille_pop
	perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 4) 
	!print*, pop(i,:), perfo(i, 1)
	end do
	
	!classement de la population
	perfo_copie = perfo
	pop = classement_population(pop, perfo_copie(:,1))
	perfo_copie = perfo
	perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)

	!Print de l'état actuel
	if (modulo(generation, 5000) == 0) then
	print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
	end if

	!génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
	if(nbre_intrus > 0) then 
		do i = taille_pop-nbre_intrus+1, taille_pop
			do j = 1, long_code
				call random(r)
				pop(i,j) = r*0.9999*18
			end do
		end do
	end if

	!Mariages
	pop = pop_mariee (pop, phi, pc, perfo)

	!mutations
	pop = pop_mutee (pop, pm)


	end do  !fin boucle sur les generations

	
	!Etait-on déjà dans le 2-gen groupe ?
	if(perfo(1,2)==1 .and. is2gen(cube_scrambled)==1) print*, 'Cube was already in the 2-gen group.'
	!A-t-on vraiment fini cette étape ?

	if( sum(perfo(1,:)) == 1700) then 
		
		print*, 'Getting into 2-gen phase is done :)'

		allocate(sequence_notation(3*perfo(1,2)))
		call ToString(pop(1,1:perfo(1,2)), sequence_notation)

		!Mise à jour de l'état de cube_scrambled si on n'y était pas déjà
		if(perfo(1,2) > 1) then
			cube_scrambled = do_sequence(cube_scrambled, pop(1,1:perfo(1,2)), perfo(1,2))
			moves_total = moves_total + perfo(1,2) !comptage pour #moves solution finale
			print*, "optimal :"
			print*, sequence_notation
			print*, "# moves"
			print*, perfo(1,2)
		end if	
	else
		print*, 'Adding 10''000 gen. to solve this part'
		GOTO 201
	end if

	print*,''

	!Prépare le print final dans le fichier
	allocate(sol2(size(sequence_notation)))
	sol2 = sequence_notation
	deallocate(pop, perfo, perfo_copie, sequence_notation)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving 2-gen !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if(twogen_bourrin > 0) then 
	print*, "calling bourrin"
	call solve_twogen_bourrin(cube_scrambled, cube_solved)
	GOTO 9999
end if

300	taille_pop = 51	    !# individus
	nbre_intrus = 5		!#individus aléatoires remplacent les plus mauvais de pop à chaque gén.
	long_code  = 30		!nombre de gènes	
	Ttot = 200000 		!nombre de générations
	pm = 0.05			!proba mutation
	pc = 0.8			!proba croisement
	phi = 1.5			!poids élite
	allocate(pop(taille_pop, long_code))
	allocate(perfo(taille_pop, 2))
	allocate(perfo_copie(taille_pop,2))

	two_gen = 1

	!génération population initiale
	do i = 1, taille_pop
		do j = 1, long_code
			call random(r)
			pop(i,j) = r*0.9999*18
		end do
	end do
	

301	do generation = 1, Ttot
	
	if(two_gen) pop = mod(pop,2) + 6*(pop/6)		!2-gen-> séquence contient uniquement (F,U)
	!Raccourci des mouvements qui s'annulent
	do i = 1, taille_pop
		pop(i,:) = trim_sequence(pop(i,:))
	end do
	
	!Evaluation de la fonction objectif sur la population
	!pour fobj : 3=2x2x3, 4=is2gen, 2=entropy
	do i = 1, taille_pop
	perfo(i, :) = fobj(cube_scrambled, cube_solved, pop(i, :), long_code, 2) 
	!print*, pop(i,:), perfo(i, 1)
	end do
	
	!classement de la population
	perfo_copie = perfo
	pop = classement_population(pop, perfo_copie(:,1))
	perfo_copie = perfo
	perfo = classement_population(perfo, perfo_copie(:,1))  !subtil :)

	!Print de l'état actuel
	if (modulo(generation, 10000) == 0) then
	print*, "Generation / elite/ moy. ", generation, perfo(1,:), float(sum(perfo(:,1))/taille_pop)
	end if

	!génocide des nmbre_intrus plus mauvais de la pop et remplacement par des aléatoires
	if(nbre_intrus > 0) then 
		do i = taille_pop-nbre_intrus+1, taille_pop
			do j = 1, long_code
				call random(r)
				pop(i,j) = r*0.9999*18
			end do
		end do
	end if

	!Mariages
	pop = pop_mariee (pop, phi, pc, perfo)

	!mutations
	pop = pop_mutee (pop, pm)


	end do  !fin boucle sur les generations

	

	!Si on a fini, on a fini, sinon on continue ... quel poète ce cyril !

	if( (perfo(1,1)+perfo(1,2))/10. == 48) then
		print*, "optimal :"
		allocate(sequence_notation(3*perfo(1,2)))
		call ToString(pop(1,1:perfo(1,2)), sequence_notation)
		print*, sequence_notation
		
		print*, "fobj / # moves / entropy"
		print*, perfo(1,:), sum(perfo(1,:))/10
		moves_total = moves_total + perfo(1,2)
		
		print*, 'Whole cube was solved in ', moves_total, ' moves :)'
	else
		print*, 'Adding 200 000 more generations to solve this stage.'
		GOTO 301
	end if
	print*,''

	!Prépare le print dans le fichier de résultat
	allocate(sol3(size(sequence_notation)))
	sol3 = sequence_notation
	deallocate(pop, perfo, perfo_copie, sequence_notation)





	!Print des résultats
999	write(10, *) 'Solution : '
	write(10, *), sol1
	write(10, *), sol2
	write(10, *), sol3
	write(10, *), 'total moves # : ', moves_total
	close(unit=10)

	Deallocate(scramble_sequence, scramble_notation, sol1, sol2, sol3)
	
9999	print*, 'End of the program.'
		read*,i
	end program RUBIK