cmaes_handlebounds.f90 8.65 KB
Newer Older
ofgeorg's avatar
ofgeorg committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
	  !-------------------------------------------------------------------------
      !  Subroutine   :          cmaes_handlebounds
      !-------------------------------------------------------------------------
      !
      !  Purpose      : This Subroutine implements boundary handling as in
      !					cmaes.m lines 728-795, it 'penalizes' function fitness
      !					with an additional term if sample was out of bounds.
      !
      !  Input		  : m				(I)	rows of arx/arxvalid
      !					n				(I)	columns of arx/arxvalid
      !					arxvalid(:,:)	(R) matrix of valid sample points
      !					arx(:,:)		(R) matrix of all sample points
      !                
      !  Remarks      : The technical details of boundary handling are not
      !					completely understood yet, so this is just a 1:1
      !					translation of Matlab Code (cmaes.m lines 728ff).
      !					bnd.flgscale as it appears in Matlab code, is not imple-
      !					mented as it is once set to 0 and never changed afterw.
      !
      !  References   : 
      !					@article{Hansen:2007,
	  !					Author = {Hansen, Nikolaus},
	  !					Keywords = {CMA},
	  !					Title = {{The CMA Evolution Strategy: A Tutorial}},
	  !					Year = {2007}}
	  !					
      !  Revisions    :
      !-------------------------------------------------------------------------
      !-------------------------------------------------------------------------
30 31 32 33
      !  pCMALib: a parallel fortran 90 library for the evolution strategy with
      !           covariance matrix adaptation
      !  Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
      !  MOSAIC group, ETH Zurich, Switzerland
ofgeorg's avatar
ofgeorg committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
      !-------------------------------------------------------------------------
      SUBROUTINE cmaes_handlebounds(arxvalid,arx,m,n)
      !-------------------------------------------------------------------------
	  !  Modules
	  !-------------------------------------------------------------------------
	  USE cmaes_param_mod
      USE cmaes_mod
      USE cmaes_opts_mod
      USE tool_mrgrnk_mod
      IMPLICIT NONE
      !-------------------------------------------------------------------------
	  !  Parameters
	  !-------------------------------------------------------------------------      
      INTEGER,INTENT(in)								:: m,n
      REAL(MK),DIMENSION(m,n),INTENT(in)				:: arxvalid
      REAL(MK),DIMENSION(m,n),INTENT(in)				:: arx
      !-------------------------------------------------------------------------
	  !  Local Variables
	  !-------------------------------------------------------------------------      
      REAL(MK),DIMENSION(2)								:: val
      REAL(MK)											:: value ! temp
      REAL(MK),DIMENSION(input%N)						:: diag
      REAL(MK)											:: mean
      INTEGER											:: i,k,maxI
      LOGICAL,DIMENSION(size(bnd%dfithist))				:: mask
      LOGICAL,DIMENSION(input%N)						:: idx
      REAL(MK),DIMENSION(input%N)						:: tx
      INTEGER,DIMENSION(size(bnd%dfithist))				:: dfitidx
      REAL(MK),DIMENSION(size(bnd%dfithist))			:: dfitsort
      REAL(MK),DIMENSION(2)								:: prct
      
	  dfitsort = 0.0_MK
	  dfitidx = 0
	  prct = (/ 25.0_MK,75.0_MK /)
	  maxI = 0

      !-------------------------------------------------------------------------
	  !  Initialize matrix mask
	  !-------------------------------------------------------------------------
      DO i = 1, size(bnd%dfithist)
      	IF(bnd%dfithist(i) .NE. 0.) THEN
      		mask(i) = .TRUE.
      		maxI = maxI + 1 ! sum up non-zero elements of bnd%dfithist vector
      	ELSE
      		mask(i) = .FALSE.
      	END IF
      END DO
      
      !-------------------------------------------------------------------------
	  !  Get main diagonal of Covariance matrix C
	  !-------------------------------------------------------------------------
	  DO i = 1, input%N
	  	diag(i) = C(i,i)
	  END DO

	  !-------------------------------------------------------------------------
	  !  ...and calculate mean value
	  !-------------------------------------------------------------------------
      mean = sum(diag)/real(input%N)

      !-------------------------------------------------------------------------
  	  !  Get delta fitness values
  	  !-------------------------------------------------------------------------
  	  CALL cmaes_myprctile(fitness%raw, lambda, prct, 2, val)
  	  value = (val(2) - val(1)) / input%N / mean /sigma**2
  	  
  	  !-------------------------------------------------------------------------
  	  !  Catch non-sensible values
  	  !-------------------------------------------------------------------------
  	  IF(value .GE. posInf) THEN
#ifdef __HAVE_MPI__
	  WRITE(*,*) 'Process ', MY_RANK, ' warning: Non-finite fitness range'
#else		  	  
  	  	WRITE(*,*) 'Warning: Non-finite fitness range'
#endif  	
  		value = maxval(bnd%dfithist)
  	  ELSE IF(value .EQ. 0.) THEN	! Happens if all points are out of bounds
  		value = minval(bnd%dfithist,mask)
  	  ELSE IF(bnd%validfitval .EQ. 0) THEN	! First sensible val
  		bnd%dfithist = 0.
  		bnd%validfitval = 1
      END IF
      
      !-------------------------------------------------------------------------
  	  !  Store delta fitness values
  	  !-------------------------------------------------------------------------
      bndd:DO i = 1, size(bnd%dfithist)
      	IF(bnd%dfithist(i) .EQ. 0) THEN	! Store at first unassigned position
      		bnd%dfithist(i) = value
      		EXIT bndd
      	ELSE IF(i .EQ. size(bnd%dfithist)) THEN		! If all positions are
      		DO k = 1, (size(bnd%dfithist)-1)		! assigned, then shift
      			bnd%dfithist(k) = bnd%dfithist(k+1) ! values to the left
      		END DO					! and store new value at last position
      		bnd%dfithist(i) = value
      	END IF
      END DO bndd
      CALL cmaes_xintobounds(options%LBounds,options%Ubounds,input%N,xmean,tx,idx)
      
      !-------------------------------------------------------------------------
  	  !  Set initial weights
  	  !-------------------------------------------------------------------------
  	  IF(bnd%iniphase) THEN
  	  	IF(any(idx)) THEN
  	  	  ! Set i to the number of non-zero elements of the bnd%dfithist vector	
  	  	  i = maxI
  	  	  IF(i .EQ. 1) THEN
  	  	  	value = bnd%dfithist(1)
  	  	  ELSE								! Median is saved in 'value'
			!sort bnd%dfithist in ascending order
			CALL mrgrnk(bnd%dfithist(1:i),dfitidx)
			DO k = 1, i
				dfitsort(k) = bnd%dfithist(dfitidx(k))
			END DO
			!and get median value
			IF(mod(i,2) .EQ. 0) THEN
				value = (dfitsort(i/2)+dfitsort(i/2+1))/2.0_MK
			ELSE
				value = dfitsort((i-1)/2+1)
			END IF  	  	
	      END IF
	      
	      WHERE(bnd%isbounded)
	      	diag = diag/mean
	      	bnd%weights = 2.0002_MK * value / diag
	      END WHERE
	      IF((bnd%validfitval+countIter) .GT. 2) bnd%iniphase = .FALSE.
	    END IF  
      END IF
      
      !-------------------------------------------------------------------------
  	  !  Increase weights
  	  !-------------------------------------------------------------------------
      IF(any(idx)) THEN
      	!-----------------------------------------------------------------------
  	    !  Judge distance of xmean to boundary
  	    !-----------------------------------------------------------------------
      	tx = xmean - tx
      	DO i = 1, input%N
      		idx(i) = ((idx(i)) .AND. (abs(tx(i)) .GT. 3.0_MK* &
      			max(1.0_MK,sqrt(real(input%N))/mueff) * sigma * sqrt(diag(i))))
      		idx(i) = ( idx(i) .AND. &
      			(sign(1.0_MK,tx(i)) .EQ. sign(1.0_MK,(xmean(i)-xold(i)))) )
      	END DO
      	!-----------------------------------------------------------------------
  	    !  Only increase if xmean is moving away
  	    !-----------------------------------------------------------------------
  	    WHERE(idx)
  	    	bnd%weights = 1.2_MK**(max(1.0_MK,mueff/10.0_MK/real(input%N))) &
  	    		* bnd%weights
  	    END WHERE
      END IF
      
      !-------------------------------------------------------------------------
  	  !  Calculate scaling biased to unity, product is one
  	  !-------------------------------------------------------------------------
  	  !mean = sum(log(diag))/real(input%N)
      !bnd%scales = exp(0.9*(log(diag)) - mean)

      !-------------------------------------------------------------------------
  	  !  Assigned penalized fitness
  	  !-------------------------------------------------------------------------
      bnd%arpenalty = matmul( (bnd%weights / bnd%scales), (arxvalid - arx)**2 )
	  fitness%sel = fitness%raw + bnd%arpenalty   

	  RETURN
	  END SUBROUTINE cmaes_handlebounds