C--   A single instance of OPTICOOL's BOSS program controls the entire
C--   optimization process, and utilizes the efforts of many Worker
C--   processors to perform the required simulations using ICOOL or
C--   similar optimization packages.

C--   BOSS 13.0  August 16, 2001

      program Boss

      implicit none

C--   Include the filename and unit number definitions
      include 'boss.inc'



C--   Maximum number of simplex vertices
      integer maxVerts
      parameter (maxVerts = maxOptVars + 1)

C--   Actual number of simplex vertices, always one more than the
C--   number of optimization parameters
      integer nVerts

C--   Initial simplex vertex coordinate
      real initSimpBase
      parameter (initSimpBase = 0.1)

C--   Initial simplex size
      real initSimpSize
      parameter (initSimpSize = 0.1)


C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

C--   Simplex vertices, by coordinate index and vertex index
C--   If there are N optimization parameters, then each vertex has
C--   N coordinates, and there are N+1 vertices in each simplex.
      real vertex (maxOptVars, maxVerts)

C--   Indexes over coordinate number (ic) and vertex number (iv)
      integer ic, iv

C--   Indexes over pass number
      integer ip, ipc

C--   Coordinates of pBar, the centroid of the non-worst vertices
      real pBar (maxOptVars)

C--   Coordinates of the reflection vector (pBar - worstVertex)
      real rVector (maxOptVars)

C--   Coordinates of the reflection vertex
      real rVertex (maxOptVars)

C--   Coordinates of the extended or contracted vertex
      real ecVertex (maxOptVars)

C--   Retained counter for each vertex in the simplex
      integer nRetained (maxVerts)

C--   Simulated and total merits for each simplex vertex, by vertex number
C--   In the real world, simulator merit is derived from simulation merit
C--   variables, whereas the total merit also includes pseudo-merits
C--   derived from parameter boundary conditions. The optimization
C--   algorithm uses the total merit, but also reports the simulated
C--   merit.
      real meritSim (maxVerts)
      real merit (maxVerts)

C--   Vertex index of the best, worst, 2nd-worst vertices in the simplex
      integer ivBest, ivWorst, iv2nd

C--   Merits of the best, worst, 2nd worst vertices
      real meritBest, meritWorst, merit2nd

C--   Merits of the reflection vertex
      real meritSimRefl, meritRefl

C--   Merits of the extension/contraction vertex
      real meritSimExtn, meritExtn

C--   Optimization pass number
      integer passCount

C--   Character string to hold the strategy descriptor
      character*60 strategy

C--   Center of gravity of the simplex for each pass
      real centGrav (maxOptVars, maxPass)

C--   Five-pass average of simplex centers of gravity
      real centGravAvg (maxOptVars)

C--   Center-of-gravity distances to five-pass average
      real cgDist (5)

C--   Number of CGs to far from the five-pass average
      integer nTooFar

C--   Stop flag -- true if this should be the last optimization pass.
      logical stopFlag


C--   Declare external function types
      real MeritFunc
      real UserOptValue

C--   Worker number
      integer iWorker

C--   Timekeeping variables
      character*30  CTime, now
      integer*8 Time8
      real runTime, Second

C--   Index over summary configurations
      integer isc

C--   Number of summary configurations for one optimization variable
      integer nsc

C--   Current value of the summary configuration's free variable
      real freeVar

C--   Spacing between adjacent summary configurations
      real spacing

C--   Number of summary simulations required
      integer nSumSim

C--   Number of summary simulations to perform this round
      integer nThisRound

C--   Value of optimization variables for a summary simulation
      real sumOptVar (maxOptVars)

C--   Index over merit variables
      integer imv


C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

C--   Zero the pass counter
      passCount = 0



C--   Open the report file; if it exists, it will be overwritten.
      open (
     &  unit=repoUnit,
     &  file=repoFile,
     &  status='unknown')

C--   Open the diagnostic file; if it exists, it will be overwritten.
      open (
     &  unit=diagUnit,
     &  file=diagFile,
     &  status='unknown')


C--   Open, read and close the optimization control file
      call ReadOptCont

C--   Open the plot file; if it exists, it will be overwritten.
      open (unit = plotUnit,
     &  file=plotFile,
     &  status='unknown')

C--   Prepare the boss-to-worker protocol. This must be done after the
C--   OptCont file has been read.
      call OpenB2W


C--   Set up the random number seeds for each Worker.
      do iWorker = 1, nWorkers
         workerSeedBase (iWorker) = -100 * iWorker
      end do

C@@   For now, the current (working) seeds are equal to the base
      do iWorker = 1, nWorkers
         workerSeedCur (iWorker) = workerSeedBase (iWorker)
      end do



C--   Number of vertices = number of optimization variables + 1
      nVerts = nOptVars + 1

C--   Write the fixed information; this will change.
      write (repoUnit,800) nOptVars, nVerts, initSimpBase, initSimpSize
      write (diagUnit,800) nOptVars, nVerts, initSimpBase, initSimpSize
 800  format (/'Variable-size Simplex Optimizer'/
     &         ' Number of parameters . . . . . ', i6/
     &         ' Number of vertices in simplex  ', i6/
     &         ' Initial simplex base . . . . . ', f6.2/
     &         ' Initial simplex size . . . . . ', f6.2/)

      call TimeOut ('Starting optimization', -9999)

C--   Set the initial values of the first nVerts vertices

C--   The first (base) vertex is at (initSimpBase, initSimpBase . . .)
      do  ic = 1, nOptVars
         vertex (ic, 1) = initSimpBase
      end do
      nRetained (1) = 0

C--   The remaining vertices are a distance initSimpSize from the
C--   base vertex; simplex edges are parallel to axis directions
      do  iv = 2, nVerts
         do  ic = 1, nOptVars
            if (ic .eq. iv-1) then
               vertex (ic, iv) = vertex (ic, 1) + initSimpSize
            else
               vertex (ic, iv) = vertex (ic, 1)
            end if
         end do
         nRetained (iv) = 0
      end do

C##   The initial simplex vertices have now been defined. Next, we must
C##   compute the merit for all of the vertices of the simplex.
      do  iv = 1, nVerts
         meritSim (iv) = MeritFunc (nOptVars, vertex(1,iv))
         merit (iv)    = meritSim (iv)
      end do


C--   Compare the merits of the various vertices . . .
C--   The vertex with the best merit
      ivBest = 1
      meritBest = merit (ivBest)

      do  iv = 2, nVerts
         if (merit (iv) .gt. meritBest) then
            ivBest = iv
            meritBest = merit (ivBest)
         end if
      end do


C--   The vertex with the worst merit
      ivWorst = 1
      meritWorst = merit (ivWorst)

      do  iv = 2, nVerts
         if (merit (iv) .lt. meritWorst) then
            ivWorst = iv
            meritWorst = merit (ivWorst)
         end if
      end do

C--   The vertex with the second-worst merit
      iv2nd = 1
      if (iv2nd .eq. ivWorst) iv2nd = 2
      merit2nd = merit (iv2nd)

      do  iv = 1, nVerts
         if (iv .ne. ivWorst) then
            if (merit (iv) .lt. merit2nd) then
               iv2nd = iv
               merit2nd = merit (iv2nd)
            end if
         end if
      end do



C---- OPTIMIZATION PASS LOOP
C--   Before coming here you must have:
C--     A completely filled in vertex coordinate table for a simplex
C--     Merits computed for each vertex
C--     Valid retained counters for each vertex
C--     The best, 2nd-worst and worst vertices defined

 100  continue
      passCount = passCount + 1

C--   Generate a new simplex. First, find pBar(), the centroid
C--   (component-by-component average) of all vertices except the worst.

      do  ic = 1, nOptVars
         pBar (ic) = 0.
         do iv = 1, nVerts
            if (iv .ne. ivWorst) then
               pBar (ic) = pBar (ic) + vertex (ic, iv)
            end if
         end do
         pBar (ic) = pBar (ic) / (nVerts - 1)
      end do

C--   Compute the reflection vector, tail at the worst vertex, head at
C--   the centroid. rVector = pBar - Worst

      do  ic = 1, nOptVars
         rVector (ic) = pBar (ic) - vertex (ic, ivWorst)
      end do

C--   Compute the default reflection vertex -- the prime candidate for
C--   being the new vertex in the next simplex. rVertex = pBar + rVector.

      do  ic = 1, nOptVars
         rVertex(ic) = pBar (ic) + rVector (ic)
      end do

C##   Compute the merit of the reflection vertex
      meritSimRefl = MeritFunc (nOptVars, rVertex)
      meritRefl    = meritSimRefl

C--   If the merit at the reflection vertex is between the 2nd-worst
C--   and the best vertices in the current simplex, the new simplex
C--   is generated from the old one by throwing out the worst vertex
C--   and inserting the reflection vertex.

      if (meritRefl .ge. merit (iv2nd) .and.
     &    meritRefl .le. merit (ivBest))  then

         do  iv = 1, nVerts
            nRetained (iv) = nRetained (iv) + 1
         end do

         do  ic = 1, nOptVars
            vertex (ic, ivWorst) = rVertex (ic)
         end do

         nRetained (ivWorst) = 0

         meritSim (ivWorst) = meritSimRefl
         merit (ivWorst) = meritRefl
         strategy = 'Reflection -- no additional vertices tested'
       
C--   If the merit at the reflection vertex is very high -- even better
C--   that the best vertex in the current simplex, try going further
C--   in the same direction.  ecVertex = rVertex + rVector

      else  if (meritRefl .gt. merit (ivBest))  then

         do  ic = 1, nOptVars
            ecVertex (ic) = rVertex (ic) + rVector (ic)
         end do

C##      Compute the merit of the extension vertex
         meritSimExtn = MeritFunc (nOptVars, ecVertex)
         meritExtn    = meritSimExtn

C--      If the extension vertex is better than the old best . . .
         if (meritExtn .ge. merit (ivBest)) then

C--         In this case we go with the extension vertex. The new
C--         simplex will be made from the old one by throwing out the
C--         worst vertex and adding the extension vertex.

            do  iv = 1, nVerts
               nRetained (iv) = nRetained (iv) + 1
            end do

            do  ic = 1, nOptVars
               vertex (ic, ivWorst) = ecVertex (ic)
            end do

            nRetained (ivWorst) = 0

            meritSim (ivWorst) = meritSimExtn
            merit (ivWorst) = meritExtn
            strategy = 'Extension -- simplex stretched'

C--      However, if the extension vertex falls short of the old best . . .
         else

C--         In this case, we stay with the standard reflection vertex.
C--         As before, make up the new simplex from the old one by
C--         replacing the worst vertex with the reflection vertex.

            do  iv = 1, nVerts
               nRetained (iv) = nRetained (iv) + 1
            end do

            do  ic = 1, nOptVars
               vertex (ic, ivWorst) = rVertex (ic)
            end do

            nRetained (ivWorst) = 0

            meritSim (ivWorst) = meritSimRefl
            merit (ivWorst) = meritRefl
            strategy = 'Reflection after extension attempt failed'

         end if


C--   Finally, we must consider the cases in which the reflection vertex
C--   is pretty crummy -- worse than the 2nd worst vertex in the existing
C--   simplex, perhaps even worse than the very worst. In this case we
C--   will contract the simplex.

      else

C--      If the reflection vertex isn't worse than the worst . . .
         if (meritRefl .ge. merit (ivWorst)) then

C--         Compute the coordinates for the contraction vertex in
C--         the forward direction
            do  ic = 1, nOptVars
               ecVertex (ic) = pBar (ic) + 0.5 * rVector (ic)
            end do

            strategy = 'Contraction in the forward direction'

C--      If the reflection vertex is the very worst . . .
         else

C--         Compute the coordinates for the contraction vertex in
C--         the backward direction
            do  ic = 1, nOptVars
               ecVertex (ic) = pBar (ic) - 0.5 * rVector (ic)
            end do

            strategy = 'Contraction in the backward direction'

         end if

C--      Stuff the contraction vertex into the new simplex; calculate
C--      the merit and stuff it in too. Don't forget the retained counters.
C--      The new simplex is made from the old one by throwing out the
C--      worst vertex and replacing it with the (forward or backward)
C--      contracted vertex.

         do  iv = 1, nVerts
            nRetained (iv) = nRetained (iv) + 1
         end do

         do  ic = 1, nOptVars
            vertex (ic, ivWorst) = ecVertex (ic)
         end do

         nRetained (ivWorst) = 0

C##      Compute the merit of the contraction vertex
         meritSimExtn = MeritFunc (nOptVars, ecVertex)
         meritExtn    = meritSimExtn

         meritSim (ivWorst) = meritSimExtn
         merit (ivWorst) = meritExtn

      end if

C--   At this point we have a brand new simplex. The coordinates are
C--   entered in the vertex coordinate table. The merits for all of
C--   vertices are up to date. The retained counts have been updated.
C--   However, note that pointers like ivWorst and ivBest are no longer
C--   valid.


C--   Compute the center of gravity of the simplex for this pass
      do  ic = 1, nOptVars
      centGrav (ic, passCount) = 0.
         do  iv = 1, nVerts
            centGrav (ic, passCount) =
     &       centGrav (ic, passCount) + vertex (ic, iv)
         end do
         centGrav (ic, passCount) = centGrav (ic, passCount) / nVerts
      end do

C--   If we've done at least five passes . . .
      if (passCount .ge. 5)  then

C--      Compute the average center of gravity for the last five passes
         do  ic = 1, nOptVars

            centGravAvg (ic) = 0.
            do  ip = passCount-4, passCount
               centGravAvg (ic) = centGravAvg (ic) + centGrav (ic, ip)
            end do
            centGravAvg (ic) = centGravAvg (ic) / 5.

         end do

C--      Compute the distance of each center of gravity from the average
         ipc = passCount-4
         do  ip = 1,5
            cgDist (ip) = 0.
            do  ic = 1, nOptVars
               cgDist (ip) = cgDist (ip) +
     &          (centGravAvg (ic) - centGrav (ic, ipc))**2
            end do
            cgDist (ip) = sqrt (cgDist (ip))
            ipc = ipc + 1
         end do

C--      If all of the last five centers of gravity lie within
C--      stopCriterion of the five-pass average, set the stop flag
         nTooFar = 0
         do  ip = 1, 5
            if (cgDist (ip) .gt. stopCriterion)  nTooFar = nTooFar + 1
         end do

         stopFlag = nTooFar .lt. 1

C--   We haven't done 5 passes, so it's too early to decide to stop
      else

         stopFlag = .false.

      end if

      write (repoUnit,850) passCount, stopFlag
      write (diagUnit,850) passCount, stopFlag
      write (*,850)        passCount, stopFlag
 850  format (/'Report for optimizer pass ', i3,
     &         '   Stop Flag = ', l1)

      write (repoUnit,851) strategy
      write (diagUnit,851) strategy
      write (*,851)        strategy
 851  format (/'Strategy:  ', a)

      do  iv = 1, nVerts
         write (repoUnit,852)
     &    merit (iv), (vertex(ic,iv), ic=1,nOptVars)
         write (diagUnit,852)
     &    merit (iv), (vertex(ic,iv), ic=1,nOptVars)
         write (*,852) merit (iv), (vertex(ic,iv), ic=1,nOptVars)
 852     format (f6.3, ' | ', 8f8.3)
      end do

      write (repoUnit,853) (centGrav (ic,passCount), ic=1,nOptVars)
      write (diagUnit,853) (centGrav (ic,passCount), ic=1,nOptVars)
      write (*,853) (centGrav (ic,passCount), ic=1,nOptVars)
 853  format ('centGrav ', 8f8.3)


C--   Stop after the maximum pass count is exceeded or the stop
C--   flag is set.
      if ((passCount .ge. maxOptSims) .or. stopFlag)  go to 200



C--   Characterize the vertices that might be part of the next simplex.
C--   The 2nd worst vertex from the previous simplex is ALWAYS the
C--   "worst" (next-to-be-discarded) vertex from the last one. This
C--   is a strange convention explained in Walters et al.

      ivWorst = iv2nd

C--   The vertex with the second-worst merit
      iv2nd = 1
      if (iv2nd .eq. ivWorst) iv2nd = 2
      merit2nd = merit (iv2nd)

      do  iv = 1, nVerts
         if (iv .ne. ivWorst) then
            if (merit (iv) .lt. merit2nd) then
               iv2nd = iv
               merit2nd = merit (iv2nd)
            end if
         end if
      end do

C--   The vertex with the best merit
      ivBest = 1
      meritBest = merit (ivBest)

      do  iv = 2, nVerts
         if (merit (iv) .gt. meritBest) then
            ivBest = iv
            meritBest = merit (ivBest)
         end if
      end do

C--   Next optimization pass . . .
      go to 100



C--   CHARACTERIZE THE OPTIMUM . . .
 200  continue

C--   Store the last center of gravity as the optimum
      do  ic = 1, nOptVars
         OptimumNorm (ic) = centGrav (ic, passCount)
      end do

C--   Convert the optimum values to user coordinates
      do  ic = 1, nOptVars
         OptimumUser (ic) = UserOptValue (ic, OptimumNorm(ic))
      end do

C--   Type out the location of the optimum
      write (repoUnit, 860)
      write (diagUnit, 860)
      write (*, 860)
 860  format (//'LOCATION OF THE OPTIMUM - - -'/
     &          '  Optimization variable name'/
     &          '  Optimal value'/
     &          '  Normalized optimal value'/
     &          '  Lower bound'/
     &          '  Upper bound')

      do  ic = 1, nOptVars

         write (repoUnit, 861) optVarName(ic)(1:optVarNameLong),
     &    OptimumUser (ic), OptimumNorm (ic), OptVarLo (ic),
     &    OptVarHi (ic)
         write (diagUnit, 861) optVarName(ic)(1:optVarNameLong),
     &    OptimumUser (ic), OptimumNorm (ic), OptVarLo (ic),
     &    OptVarHi (ic)
         write (*, 861) optVarName(ic)(1:optVarNameLong),
     &    OptimumUser (ic), OptimumNorm (ic), OptVarLo (ic),
     &    OptVarHi (ic)
 861     format (/a, 4f12.5)

      end do

C--   Using information from OptCont, prepare a list of all
C--   configurations to be simulated for the summary.

C--   Zero the summary configuration counter
      nSumConfigs = 0

C--   For each optimization variable . . .
      do  ic = 1, nOptVars

C--      Number of summary configurations for this variable
         nsc = optVarSumConfigs (ic)

C--      Spacing between configurations (evenly spaced in (0,1))
         spacing = 1. / (float(nsc) - 1.)

C--      Zero the value of the free variable
         freeVar = 0.

C--      For each new summary configuration . . .
         do  isc = 1, nsc

C--         Don't let the table overflow
            if (nSumConfigs .ge. maxSumConfigs) STOP

C--         Increment the summary configuration number
            nSumConfigs = nSumConfigs + 1

C--         Record the index of the free variable
            sumFreeIndex (nSumConfigs) = ic

C--         Make the last configuration for each variable exactly 1
            if (isc .eq. nsc)  freeVar = 1.0

C--         Record the value of the free variable
            sumFreeValue (nSumConfigs) = freeVar

C--         Next value of the free variable
            freeVar = freeVar + spacing

         end do

      end do

C--   Output the summary configuration list
      write (diagUnit, 870)
 870  format (//'Summary configurations to be computed - - -'/)
      do  isc = 1, nSumConfigs
         write (diagUnit, 871) isc, sumFreeIndex (isc),
     &    sumfreeValue (isc)
 871     format (5x, i3, 2x, i2, 2x, f12.4)
      end do


C--   Initialize the counter of summary configurations simulated
      nSumSim = 0

C--   Compute the number of configurations to be simulated this
C--   round; smaller of (1) the number of workers available, and
C--   (2) the number of configuration as yet unsimulated.

      nThisRound = min (nWorkersIdle, nSumConfigs-nSumSim)


C--   Until all the summary simulations have been done . . .
      do while (nThisRound .gt. 0)

C--      For every worker . . .
         do  iWorker = 1, nWorkers

C--         If this worker is idle and ready to work . . .
            if (workerState (iWorker) .eq. 1) then

C--            Increment the summary simulation counter
               nSumSim = nSumSim + 1

C--            Prepare the optimization variables for this simulation
C--            First, load the position of the optimum
               do  ic = 1, nOptVars
                  sumOptVar (ic) = optimumNorm (ic)
               end do

C--            For the free variable, substitute the appropriate value
               sumOptVar (sumFreeIndex (nSumSim)) =
     &          sumFreeValue (nSumSim)

C--            Specify the number of particles to be propagated
               workerNPart (iWorker) = nPartSum

C--            Start the simulation; the started workers are set
C--            to the "active" state.
               call StartSim (iWorker, 0, sumOptVar)

               write (*,875) nSumSim, iWorker, SumFreeIndex (nSumSim),
     &          sumFreeValue (nSumSim)

 875           format (/'Summary Simulation - - '/
     &                  '  Simulation number . . . . . ', i10/
     &                  '  Worker number . . . . . . . ', i10/
     &                  '  Index of free variable  . . ', i10/
     &                  '  Value of free variable  . . ', f10.4)

C--            Record the worker task so we know where to store the
C--            results of the simulation
               workerTask (iWorker) = nSumSim

            end if

C--         If all the needed simulations are underway, leave the
C--         remaining idle workers to stay idle.
            if (nSumSim .ge. nSumConfigs) go to 300

         end do


C--      All needed simulation for this round have been started
 300     continue

C--      While there are any active workers . . .
         do while (nWorkersActive .gt. 0)

C--         Wait until some active worker reports results
            call WaitResults (iWorker)

C--         Fetch the results from the worker's Merit Out file
C--         Results are stored in the meritVarValue array.
            call FetchResults (iWorker)

C--         Compute the merits of this worker's configuration
            call ComputeWorkerMerit (iWorker)

C--         Move merit-related data into the summary arrays
            isc = workerTask (iWorker)
            do  imv = 1, nMeritVars
               sumMeritVar   (isc, imv) = meritVarValue (imv, iWorker)
               sumMeritComp  (isc, imv) = meritContrib  (imv, iWorker)
            end do

            sumMeritTotal (isc) = meritWorker (iWorker)

C--         Put the worker back in the idle state
            workerState (iWorker) = 1
            nWorkersDone = nWorkersDone - 1
            nWorkersIdle = nWorkersIdle + 1

         end do

C--      Number of simulations to perform next round
         nThisRound = min (nWorkersIdle, nSumConfigs-nSumSim)

      end do


C--   All summary simulations completed; prepare the final reports

      write (*, 879)
 879  format (//'Summary simulations completed'/)

      write (diagUnit, 880) nSumConfigs
 880  format (///'Data for ', i3, ' summary configurations'/)
      do  isc = 1, nSumConfigs
         write (diagUnit, 881) isc, sumFreeIndex (isc),
     &    sumFreeValue (isc), sumMeritTotal (isc)
 881     format (/i3, ': ', i3, f12.4, f12.4)
         write (diagUnit, 882)
     &    (sumMeritVar (isc, imv), imv=1,nMeritVars)
 882     format (6f12.4)
         write (diagUnit, 883)
     &    (sumMeritComp (isc, imv), imv=1,nMeritVars)
 883     format (6f12.4)
      end do

C--   Set up the display descriptors (page, row, column).


C--   Output the summary plot postscript page descriptors.


C--   Output the summary characterizations as text tables
      call printSumConfigs (repoUnit)
      call printSumConfigs (diagUnit)
      call printSumConfigs (6)

      call printSumConfigs2 (repoUnit)
      call printSumConfigs2 (diagUnit)
      call printSumConfigs2 (6)



C--   All done; shut down and exit.
 400  continue

C--   Send "fall on your sword" messages to the Workers
      call KillWorkers

C--   Output the current date/time
      call TimeOut ('Ending Optimization', -9999)

C--   Output the CPU time used by the Boss
      runTime = Second ()
      write (*,894) runTime
      write (repoUnit,894) runTime
      write (diagUnit,894) runTime
 894  format (/'Boss CPU Time (seconds): ', f10.3)

C--   Shut down the Boss-to-Worker connection
      call CloseB2W

C--   Close the plot file
      close (unit = plotUnit)

C--   Close the report file
      close (unit = repoUnit)

C--   Close the diagnostic file
      close (unit = diagUnit)

      stop

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      real function MeritFunc (nParams, params)

      implicit  none

      integer nParams
      real params (nParams)

      include 'boss.inc'


C--   Declare function types
      real UserOptValue

C--   General purpose index variables
      integer i, j, k

C--   Diagnostic loop index
      integer ii

C--   Worker index
      integer iWorker

C--   Registers for average and std. dev. calculations
      real sumX, sumX2, nf, sdv2

C--   Boundary-related merit for the configuration
      real meritBound

C--   Contribution to the boundary-related merit by one parameter
      real bMerit

C--   Final merit calculation -- simulation and bounds
      real meritFinal

C--   Upper and lower bounds for a given parameter
      real boundLo, boundHi

C--   Normalized parameter value,  boundary merit calculations
      real paramNorm

C--   Current value of parameter(i)
      real paramValue

C--   User-defined boundary steepness parameter
      real steep




C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      write (diagUnit, 700) nParams
 700  format (//'Entering routine MeritFunc - -' /
     &          '  Number of parameters = ', i2)
      do  ii = 1, nParams
         write (diagUnit, 701) ii, params(ii)
 701     format ('  Parameter ', i2, ':  ', f10.4)
      end do
      write (diagUnit, 702)
 702  format ('Calling subroutine Simulate'/)

      write (*, 800)
      write (repoUnit, 800)
 800  format (/'Simulating configuration and assessing merit - -')
      do  i = 1, nParams
         write (*,801) i, optVarName (i)(1:optVarNameLong),
     &    UserOptValue (i, params(i)), params(i)
         write (repoUnit,801) i, optVarName (i)(1:optVarNameLong),
     &    UserOptValue (i, params(i)), params(i)
 801     format (3x, i2, '. ', a, 2x, f10.4, ' (', f7.4, ' )' )
      end do


C--   Compute the boundary merits for this configuration

C--   Start out assuming that all parameters are in-bounds
      meritBound = 1.

C--   For every optimization variable . . .
      do  i = 1, nParams

C--      Current value of this optimization variable
         paramValue = params (i)

C--      Lower bound and upper bound
         boundLo = optVarLo (i)
         boundHi = optVarHi (i)

C--      Map the value onto the interval [0,1]
C@@      paramNorm = (paramValue - boundLo) / (boundHi - boundLo)
         paramNorm = paramValue

C--      If the variable is above the upper bound . . .
         if (paramNorm .gt. 1.0001)  then

C--         User-specified steepness of the boundary merit
            steep = optVarSteep (i)

C--         Calculate the boundary merit
            bMerit = exp ( -(steep**2 * (paramNorm-1.)**2.5))

C--      If the variable is below the lower bound . . .
         else  if (paramNorm .lt. -0.0001) then

C--         User-specified steepness of the boundary merit
            steep = optVarSteep (i)

C--         Calculate the boundary merit
            bMerit = exp ( -(steep**2 * (-paramNorm)**2.5))

C--      If the variable is (almost) within bounds . . .
         else

C--         The boundary merit is 1
            bMerit = 1.0

         endif

C--      Multiply this into the overall boundary merit
         meritBound = meritBound * bMerit

      end do
c
c         write(*,902) meritBound, bMerit
c902   format('meritBound, bMerit', 2f12.4)
c
C@@   Eventually the boundary merit will be used to simulate
C@@   lightly or not at all those configurations with outrageously
C@@   low boundary merits.



C--   Perform the simulations - - -
C--    Prepare input files for the workers
C--    Tell the workers to get busy
C--    Wait for the returning merit files
C--    Store the merit variables returned from each worker


C--   Divide the workload among the workers according to speed
      call AssignWorkload (nPartSim)

C--   Start the simulations for this configuration on all idle workers
      do  iWorker = 1, nWorkers
         if (workerState (iWorker) .eq. 1) then
            call StartSim (iWorker, 0, params)
         end if
      end do


C--   Fetch the simulation results and process them
 10   do while (nWorkersActive .gt. 0)

C--      Wait for some active worker to report results
         call WaitResults (iWorker)

C--      Fetch and store the results from the worker's merit file
         call FetchResults (iWorker)

C--      Compute the configuration merit based on this worker's
C--      simulation
         call ComputeWorkerMerit (iWorker)

      end do


C--   Compute the average merit over all participating workers
      if (nWorkers .eq. 1) then
         meritAvg = meritWorker(1)
         meritSdv = 0.
      else
         sumX = 0.
         sumX2 = 0.
         do  iWorker = 1, nWorkers
            sumX = sumX + meritWorker(iWorker)
            sumX2 = sumX2 + meritWorker(iWorker)**2
         end do
         nf = float (nWorkers)
         meritAvg = sumX / nf
         sdv2 = (sumX2 - sumX**2 / nf) / (nf - 1.)
         sdv2 = max (sdv2, 0.)
         meritSdv = sqrt (sdv2)
      end if


C--   The final merit is the average merit from the simulations
C--   multiplied by the boundary merit
      meritFinal = meritBound * meritAvg

C--   The function returns the total merit, including both
C--   boundary merits and simulation-derived merits.
      MeritFunc = meritFinal


C--   Set all the participating Workers back to "idle" state
      do  iWorker = 1, nWorkers
         if (workerState (iWorker) .eq. 3) then
            workerState (iWorker) = 1
            nWorkersDone = nWorkersDone - 1
            nWorkersIdle = nWorkersIdle + 1
         end if
      end do

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine OpenB2W
      implicit none

      include 'boss.inc'

C--   Initialize the boss-to-worker protocol

      integer istat
      integer Remove
      integer iWorker

C--   Make sure there are no old results files lying around to confuse
C--   startup.
      call System ('rm -f ~/bd/MeritOut.??', istat)
      call System ('rm -f ~/bd/RESULTS.??', istat)

C--   Initialize the Worker state counters
      nWorkersIdle = nWorkers
      nWorkersActive = 0
      nWorkersDone = 0

C--   Set all Workers to "idle"
      do  iWorker = 1, nWorkers
         workerState (iWorker) = 1
      end do


      write (*,820)
 820  format (//'SETTING UP THE MULTIPROCESSOR ENVIRONMENT')

      write (*,825)
 825  format (//'  Putting down any existing Worker processes')

C--   Kill all Worker processes in all processors to be used
      do  iWorker = 1, nWorkers
c         print 888, workerName(iWorker)
c888   format(a40)
         call KillProcWorkers (workerName (iWorker))
      end do

      write (*,830)
 830  format (//'  Starting new Worker processes')

C--   Start a new Worker process in each worker processor
      do  iWorker = 1, nWorkers
         call StartProcWorker (workerName (iWorker))
      end do

      write (*,835)
 835  format (//'  Transmitting BossInfo files to all Workers')

C--   Send BossInfo files out to all the Workers
      do  iWorker = 1, nWorkers

C--      Open a temporary file
         open (
     &    unit = recycUnitA,
     &    file = pathBD // 'BossInfo',
     &    status='unknown')

C--      Write the scp vector for the Merit File to be sent to
C--      the Boss
         write (recycUnitA, 800) account, bossName, relPathBD,
     &    iWorker
 800     format (a, '@', a, ':', a, '/MeritOut.', i2.2)

C--      Write the scp vector for the RESULTS file to be sent
C--      to the Boss
         write (recycUnitA, 805) account, bossName, relPathBD,
     &    iWorker
 805     format (a, '@', a, ':', a, '/RESULTS.', i2.2)

c         write(*,901) account,bossName,relPathBD,iWorker
c901      format(3(2x,a20),i6)
C--      Close the temporary file
         close (unit=recycUnitA)

C--      Transmit this file to the appropriate Worker
c         write(*,902) iWorker
c902      format('call to TransmitFile--BossInfo',i4)
         call TransmitFile (iWorker, 'BossInfo')

C--   Next worker . . .
      end do

      write (*,840)
 840  format (//'BEGINNING OPTIMIZATION'//)

      return
      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine CloseB2W
      implicit none

      include 'boss.inc'

C--   Terminate the boss-to-worker protocol

      integer istat
      integer Remove

C--   Make sure there are no results file left lying around to confuse
C--   future users.
      call System ('rm -f ~/bd/MeritOut.??')
      call System ('rm -f ~/bd/RESULTS.??')


      return
      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM


      subroutine ReadOptCont
      implicit none

C--   Read the optimization control file and store the data in
C--   variables in the /OptCont/ common block.


C--   Include the common blocks and parameters
      include 'boss.inc'


C--   Unit number for reading the optimization control file
      integer optUnit
      parameter (optUnit=recycUnitA)

C--   Name of the optimization control file
      character*(*) optFile
      parameter (optFile = pathBD//'OptCont.dat')

C--   Maximum line length within the optimization control file
      integer maxLineLen
      parameter (maxLineLen = 200)

C--   One line of the input file
      character line*(maxLineLen)

C--   Length of the input line with trailing spaces trimmed
      integer lenTrim

C--   A single input character
      character*1 ich

C--   Index of leftmost non-blank and rightmost non-blank
      integer lmnb, rmnb

C--   General purpose index variables
      integer i, j

C--   Maximum number of characters in a section identifier
      integer maxSectIDLen
      parameter (maxSectIDLen = 30)

C--   Section identifier, up to maxSectIDLen characters
      character*(maxSectIDLen)  sectID

C--   Number of valid section IDs
      integer nValSectID
      parameter (nValSectID = 6)

C--   Table of valid section IDs. They must be no longer than
C--   maxSectIDLen characters long. There must be exactly
C--   nValSectID of them.
      character*(maxSectIDLen)  valSectID(nValSectID)
      data valSectID /  
     & 'WorkerNames' ,
     & 'OptimizationVariables' ,
     & 'MeritVariables' ,
     & 'InputSubstitutions' ,
     & 'FileTransfers' ,
     & 'Parameters' /

C--   Current section number (1, 2 . . . nValSectID), 0 = none
      integer iSect

C--   i/o status from read operations
      integer ios

C--   Function to extract a word from a string
      logical GetWord

C--   Temporary storage for merit combination code
      character*2 comCode

C--   True if a matching entry has been found
      logical match

C--   Length of input string
      integer stringLen

C--   Optimization parameter name
      character*(maxLineLen)  paramName

C--   Optimization parameter value
      real paramValue


C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

C--   No section ID read yet
      iSect = 0

C--   Zero the section length counters
      nWorkers = 0
      nOptVars = 0
      nMeritVars = 0
      nSubVars = 0
      nInputFiles = 0

C--   Open the control file; it must exist.
      open (unit=optUnit, file=optFile, status='old')


C--   Line processing loop . . .
 100  continue

C--   Read a line from the configuration file; if end-of-file, 
C--   exit from the loop
      read (optUnit, 800, end=200) line
 800  format (a)

C--   Number of characters excluding trailing spaces
      lenTrim = Len_Trim(line)

C--   If this is a blank line, ignore it; process next line
      if (lenTrim .eq. 0)  go to 100

C--   There is at least one non-blank character. If the first
C--   non-blank character is !, this is a comment, so ignore it
C--   go back to get another line. If the first non-blank
C--   character is anything but !, continue; lmnb is the index
C--   of the leftmost non-blank.
      do  lmnb = 1, lenTrim
         ich = line(lmnb:lmnb)
         if (ich .ne. ' ') then
            if (ich .eq. '!')  then
               go to 100
            else
               go to 105
            end if
         end if
      end do
 105  continue

C--   It's not a blank line and it's not a comment, so it presumably
C--   has data that we need to read. It might be a section header
C--   (first character non-blank) or it might be an ordinary data
C--   line (indented at least one space).

C--   If this is a section header . . .
      if (lmnb .eq. 1)  then

C--      Section headers consist of a section identifier. To allow for
C--      comments following the identifier, we assume that the 
C--      identifier begins in column 1 and ends at the first space
C--      or the end of trimmed line.
         rmnb = lenTrim
         do  i = 2, lenTrim
            ich = line (i:i)
            if (ich .eq. ' ') then
               rmnb = i-1
               go to 110
            end if
         end do

C--      Extract the section identifier
 110     sectID = line(lmnb:rmnb)

C--      Lookup this section ID in the list of valid section IDs.
C--      If it's found, the sectionID number (1,2 ... nValSect) is
C--      in iSect. Until a new section header is found, all additonal
C--      lines read must belong to this section.
         do  iSect = 1, nValSectID
            if (sectID .eq. valSectID (iSect)) go to 100
         end do

C--      Unrecognized section ID; fatal error
         stop


C--   If this is not a section header . . .
      else

C--      If no section header has yet been read, this is an orphan
C--      data line; fatal error
         if (iSect .eq. 0)  stop


C--      Section 1: worker names and specifications
         if (iSect .eq. 1)  then

C--         Don't let the worker table overflow
            if (nWorkers .ge. maxWorkers)  stop

C--         Increment the worker count
            nWorkers = nWorkers + 1

C--         Extract the name of the worker processor
            if (.not. GetWord (line, lmnb, workerName (nWorkers))) stop

C--         Extract the remaining parameters from the line
            read (line (lmnb+1:), *, iostat=ios)
     &       workerSpeedNom (nWorkers)
            if (ios .ne. 0)  stop
            workerSpeedCur (nWorkers) = workerSpeedNom (nWorkers)


C--      Section 2:  optimization variable declarations
         else  if (iSect .eq. 2)  then

C--         Don't let the optimization variable table overflow
            if (nOptVars .ge. maxOptVars)  stop

C--         Increment the optimization variable counter
            nOptVars = nOptVars + 1

C--         Extract the name of the optimization variable
            if (.not. GetWord (line, lmnb, optVarName (nOptVars))) stop

C--         Extract the remaining parameters from the line
            read (line (lmnb+1:), *, iostat=ios)
     &       optVarLo (nOptVars),
     &       optVarHi (nOptVars),
     &       optVarSteep (nOptVars),
     &       optVarSumConfigs (nOptVars)
            if (ios .ne. 0)  stop


C--      Section 3:  merit variable declarations
         else  if (iSect .eq. 3)  then

C--         Don't let the merit variable table overflow
            if (nMeritVars .ge. maxMeritVars)  stop

C--         Increment the merit variable counter
            nMeritVars = nMeritVars + 1

C--         Extract the merit variable name
            if (.not. GetWord (line, lmnb, meritVarName (nMeritVars)))
     &       stop

C--         Extract the merit combination code
            if (.not. GetWord (line, lmnb, comCode)) stop

C--         Valid combine flags are M or A1 . . . A9
            ich = comCode(1:1)

C--         Combination code for M (multiply is zero)
            if (ich .eq. 'm' .or. ich .eq. 'M') then
               meritVarCombin (nMeritVars) = 0

C--         Combination code for A1 is 1 . . . A9 is 9.
            else  if (ich .eq. 'a' .or. ich .eq. 'A') then
               read (comCode(2:2), *) meritVarCombin (nMeritVars)

C--         If it's not M or A-numeral, it's a fatal error
            else
               stop

            end if

C--         Extract the remaining parameters from the line
            read (line (lmnb+1:), *) meritVarPleft (nMeritVars),
     &       meritVarPright (nMeritVars), meritVarMmin (nMeritVars),
     &       meritVarMmax (nMeritVars)



C--      Section 4: input substitution declarations
         else  if (iSect .eq. 4)  then

C--         Don't let the substitution table overflow
            if (nSubVars .ge. maxSubVars) stop

C--         Increment the substitution variable counter
            nSubVars = nSubVars + 1

C--         First word is the substitution marker
            if (.not. GetWord (line, lmnb, subMarkName (nSubVars))) stop

C--         Second word is the associated optimization variable name
            if (.not. GetWord (line, lmnb, subOptName (nSubVars))) stop

C--         Third word is the format specifier
            if (.not. GetWord (line, lmnb, subFormat (nSubVars))) stop

C--         Get the multiplier and addend from the remainder
            read (line (lmnb+1:), *, iostat=ios)
     &       subMultiplier (nSubVars), subAddend (nSubVars)
            if (ios .ne. 0)  stop


C--      Section 5: file transfer declarations
         else  if (iSect .eq. 5)  then

C--         Don't let the filename table overflow
            if (nInputFiles .ge. maxInputFiles) stop

C--         Increment the input file counter
            nInputFiles = nInputFiles + 1

C--         First word is the unmodified input file name
            if (.not. GetWord (
     &       line, lmnb, inputFilename (nInputFiles))) stop

C--         Second word (which may be absent) is the modified input
C--         file name. Make sure that if the name is absent the
C--         specifier is null.
            if (.not. GetWord (
     &       line, lmnb, modifFilename (nInputFiles))) then
               modifFilename (nInputFiles) = ''
            end if

C--      Section 6: optimization parameters
         else  if (iSect .eq. 6)  then

C--         First word is the name of the parameter
            if (.not. GetWord (line, lmnb, paramName)) stop

C--         Following the name is the value
            read (line (lmnb+1:), *, iostat=ios) paramValue
            if (ios .ne. 0)  stop

C--         Search the parameter name table for this entry
            if (paramName .eq. 'nPartSim') then
               nPartSim = paramValue
            else  if (paramName .eq. 'nPartSum') then
               nPartSum = paramValue
            else  if (paramName .eq. 'stopCriterion') then
               stopCriterion = paramValue
            else  if (paramName .eq. 'maxOptSims') then
               maxOptSims = paramValue
               if (maxOptSims .gt. maxPass) maxOptSims = maxPass
            else
               STOP
            end if


C--      End of section-by-section line decoding
         end if

C--      Read another line . . .
         go to 100

      end if



C--   We've reached the end of the input file
 200  continue

C--   By now both optVars and subVars has been read, so build the index
C--   array for subVars.

C--   For every substitution declaration . . .
      do  i = 1, nSubVars

C--      Initialize search for matching optimization variable name
         j = 1
         match = .false.

C--      Try to find the cited optimization variable name
         do while (j .le. nOptVars .and. .not. match)

             if (subOptName (i) .eq. optVarName (j)) then
                subOptIndex (i) = j
                match = .true.
             else
                j = j+1
             end if

         end do

C--      If no match could be found, it's a fatal error
         if (.not. match)  stop

      end do

C--   Compute various longest-string-lengths for formatting output
      subMarkNameLong = 0
      do  i = 1, nSubVars
         stringLen = Len_Trim (subMarkName (i))
         if (stringLen .gt. subMarkNameLong)
     &     subMarkNameLong = stringLen
      end do

      subOptNameLong = 0
      do  i = 1, nSubVars
         stringLen = Len_Trim (subOptName (i))
         if (stringLen .gt. subOptNameLong)
     &     subOptNameLong = stringLen
      end do

      subFormatLong = 0
      do  i = 1, nSubVars
         stringLen = Len_Trim (subFormat (i))
         if (stringLen .gt. subFormatLong)
     &     subFormatLong = stringLen
      end do

      meritVarNameLong = 0
      do  i = 1, nMeritVars
         stringLen = Len_Trim (meritVarName (i))
         if (stringLen .gt. meritVarNameLong)
     &     meritVarNameLong = stringLen
      end do

      optVarNameLong = 0
      do  i = 1, nOptVars
         stringLen = Len_Trim (optVarName (i))
         if (stringLen .gt. optVarNameLong)
     &     optVarNameLong = stringLen
      end do

      workerNameLong = 0
      do  i = 1, nWorkers
         stringLen = Len_Trim (workerName (i))
         if (stringLen .gt. workerNameLong)
     &     workerNameLong = stringLen
      end do

      inputFilenameLong = 0
      do  i = 1, nInputFiles
         stringLen = Len_Trim (inputFilename (i))
         if (stringLen .gt. inputFilenameLong)
     &     inputFilenameLong = stringLen
      end do

      modifFilenameLong = 0
      do  i = 1, nInputFiles
         stringLen = Len_Trim (modifFilename (i))
         if (stringLen .gt. modifFilenameLong)
     &     modifFilenameLong = stringLen
      end do


C--   Close the optimization control file
      close (unit=optUnit)

C--   Print a control variable summary to the report file
      call PrintOptCont (repoUnit)
      call PrintOptCont (diagUnit)
      call PrintOptCont (6)

      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine PrintOptCont (unit)
      implicit  none

C--   Print the optimization control data to a specified unit number.
C--   The file is assumed open for writing, and is left open.



C--   Include the common blocks and parameters
      include 'boss.inc'

      integer unit

      integer i

      character*2 combinCode


C--   Worker Processors

      write (unit, 800)
 800  format (//'Worker processors - - -'/
     &          '  Internet name of the processor'/
     &          '  Nominal speed (MHz) of the processor')

      write (unit, 801) nWorkers
 801  format (/'Number of workers: ', i3/)

      do  i = 1, nWorkers
         write (unit, 805) i,
     &    workerName(i)(1:workerNameLong),
     &    workerSpeedNom (i)
      end do
 805  format (2x, i2, 2x, a, 2x, f7.1)


C--   Optimization Variables

      write (unit, 810)
 810  format (//'Optimization Variables - - -'/
     &          '  Optimization variable name'/
     &          '  Lower bound on the variable'/
     &          '  Upper bound on the variable'/
     &          '  Steepness of the boundary merit function'/
     &          '  Number of summary configurations needed')

      write (unit, 811) nOptVars
 811  format (/'Number of optimization variables: ', i3/)

      do  i = 1, nOptVars
         write (unit, 815) i,
     &    optVarName (i)(1:optVarNameLong),
     &    optVarLo (i),
     &    optVarHi (i),
     &    optVarSteep (i),
     &    optVarSumConfigs (i)
      end do
 815  format (2x, i2, 2x, a, 3f12.4, 2x, i2)


C--   Substitution Variables

      write (unit, 820)
 820  format (//'Substitution variables - - -'/
     &          '  Substitution mark name'/
     &          '  Associated optimization variable name'/
     &          '  Index of associated optimization variable'/
     &          '  Format to express substituted value with'/
     &          '  Multiplier for substitution value'/
     &          '  Addend for substitution value')

      write (unit, 821) nSubVars
 821  format (/'Number of substitution directives: ', i3/)

      do  i = 1, nSubVars
         write (unit, 825) i,
     &    subMarkName (i)(1:subMarkNameLong),
     &    subOptName (i)(1:subOptNameLong),
     &    subOptIndex (i),
     &    subFormat (i)(1:subFormatLong),
     &    subMultiplier (i),
     &    subAddend (i)
      end do
 825  format (2x, i2, 2x, a, 2x, a, ' (', i2, ')  ', a, 2f10.4)


C--   Merit Variables

      write (unit, 830)
 830  format (//'Merit variables - - -'/
     &          '  Merit variable name'/
     &          '  Combination code for this variable'/
     &          '  Variable value having 10% full merit'/
     &          '  Variable value having 90% full merit'/
     &          '  Minimum value of the merit'/
     &          '  Maximum value of the merit')

      write (unit, 831) nMeritVars
 831  format (/'Number of merit variables: ', i3/)

      do  i = 1, nMeritVars

         if (meritVarCombin (i) .eq. 0) then
            combinCode = 'M '
         else
            write (combinCode, 832) meritVarCombin (i)
 832        format ('A', i1)
         end if

         write (unit, 835) i,
     &    meritVarName (i)(1:meritVarNameLong),
     &    combinCode,
     &    meritVarPleft(i),
     &    meritVarPright(i),
     &    meritVarMmin(i),
     &    meritVarMmax(i)
 835     format (2x, i2, 2x, a, 2x, a, 4f10.4)

      end do


C--   Input file transfers

      write (unit, 840)
 840  format (//'Input files to transfer - - -'/
     &          ' Unmodified input filename'/
     &          ' Modified (substituted) input filename')

      write (unit, 841) nInputFiles
 841  format (/'Number of input files to transfer: ', i3/)

      do  i = 1, nInputFiles

         write (unit, 845) i,
     &    inputFilename (i)(1:inputFilenameLong),
     &    modifFilename (i)(1:modifFilenameLong)
 845     format (2x, i2, 2x, a, 2x, a)

      end do

C--   Optimization Parameters

      write (unit, 850) nPartSim, nPartSum, stopCriterion, maxOptSims
 850  format (//'Optimization parameters:'/
     &  '  Number of particles to propagate while optimizing: ', i10/
     &  '  Number of particles to propagate in summaries:     ', i10/
     &  '  Optimization Stop Criterion:                       ', f10.4/
     &  '  Maximum number of optimization simulations:        ', i10)


      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

C--   For each substitution marker in the specified input file,
C--   substitute an actual value (generated by the optimizer,
C--   stored in array values). Translate from optimizer coordinates
C--   to simulator coordinates. The function returns true if
C--   substitution succeeds, otherwise false.

      logical function SubstFile (inFile, outFile, values, iWorker)

      implicit none

      character*(*) inFile, outFile
      real values(*)
      integer iWorker

      include 'boss.inc'

      integer inUnit, outUnit
      parameter (inUnit = recycUnitA, outUnit = recycUnitB)

      integer ios

      character*(maxSimInLineLen) line

      logical SubstLine

c      write(*,*) 'entry SubstFile'

C--   Initialize return value to "failed"
      SubstFile = .false.

C--   open the input file
      open (unit=inUnit, file=inFile, status='old', iostat=ios)
      if (ios .ne. 0) return

C--   open the output file
      open (unit=outUnit, file=outFile, status='unknown', iostat=ios)
      if (ios .ne. 0) return


C--   Line processing loop . . .
 10   continue

C--   Read a line from the input file
      read (inUnit, 800, iostat=ios) line
 800  format (a)

C--   If it's end-of-file, close files and return success
      if (ios .eq. -1) then
         close (unit=inUnit)
         close (unit=outUnit)
         SubstFile = .true.
         return
      end if

C--   If it's an error, close files and return failure
      if (ios .ne. 0) then
         close (unit=inUnit)
         close (unit=outUnit)
         SubstFile = .false.
         return
      end if

C--   It's a good line. Try to make all substitutions required in this
C--   line. Return failure if the line substitution fails.
      if (.not. SubstLine (line, values, iWorker)) return

C--   We've received a good line with all needed substitutions made.
C--   Write it to the output file
      write (outUnit, 810, iostat=ios) line(1:Len_Trim(line))
 810  format (a)

C--   If there's any problem return failure
      if (ios .ne. 0)  return

C--   All's well; go process the next input line . . .
      go to 10

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

C--   For each substitution marker in the specified input line,
C--   substitute an actual value (generated by the optimizer,
C--   stored in array values). Translate from optimizer coordinates
C--   to simulator coordinates. The function returns true if
C--   substitution succeeds, otherwise false.

      logical function SubstLine (line, values, iWorker)

      implicit none

      character*(*) line
      real values(*)
      integer iWorker

      logical done

      logical SubstMark



C--   Initialize the return to "failure"
      SubstLine = .false.

C--   Not done yet
      done = .false.


C--   Marker processing loop . . .
      do while (.not. done)

C--      If there is at least one substitution marker in the input line,
C--      make the required substitution. If the attempt fails, return
C--      failure.
         if (.not. SubstMark (line, done, values, iWorker)) return

      end do

C--   All the markers were successfully substituted
      SubstLine = .true.
      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

C--   Find a single substitution marker in the current line..
C--   Perform the indicated substitution, returning the reconstituted
C--   line. Index points to the character following the substituted
C--   string; it is set to 0 if there are no more markers to be
C--   processed.

      logical function SubstMark (line, done, values, iWorker)

      implicit none

      character*(*) line
      logical done
      real values(*)
      integer iWorker

      include 'boss.inc'

      integer lenTrim
      integer lcb, rcb
      character*(maxSubMarkLen) mark
      integer i
      integer iopt
      real optVal
      character (maxSimInFieldLen) optValStr
      integer optValStrLen



C--   Initialize the return value to "failure"
      SubstMark = .false.

C--   Get the trimmed length of the line
      lenTrim = Len_Trim (line)

C--   Look for the leftmost {
      do  lcb = 1, lenTrim
         if  (line (lcb:lcb) .eq. '{') go to 10
      end do

C--   There are no left curly brackets. Set done and return success
      done = .true.
      SubstMark = .true.
      return

C--   We've found a left curly bracket. It must be a substitution mark.
 10   continue

C--   Find the corresponding right curly bracket
      do  rcb = lcb+1, lenTrim
         if (line (rcb:rcb) .eq. '}')  go to 20
      end do

C--   There is no corresponding right curly bracket; serious error.
      return

C--   We've found a right curly bracket too.
 20   continue

C--   Extract the substitution marker, including brackets
      mark = line (lcb:rcb)


C--   Does it stand for "number of particles to propagate"?
      if (mark .eq. '{NPART}')  then

C--      Get the number of particles to be propagated by this
C--      worker and format it appropriately
         write ( optValStr, '(i7)' ) workerNPart (iWorker)


C--   Does it stand for "seed for the random number generator"?
      else  if (mark .eq. '{RNSEED}')  then

C--      Get the random seed for this worker and format it
C--      appropriately
         write ( optValStr, '(i10)' ) workerSeedCur (iWorker)


C--   Does it stand for some substitution marker in the subMark table?
      else

C--      Look up this marker in the substitution mark table
         do  i = 1, nSubVars
            if (mark .eq. subMarkName (i)) go to 30
         end do

C--      No match; this is a bad marker. Crawl off to die.
         return

C--      This mark is one registered in the subMark tables.
 30      continue

C--      Get the associated optimization variable index
         iopt = subOptIndex (i)

C--      Get the current value of this optimization variable. The
C--      numbers in array values are normalized (in the range (0,1))
C--      and must be transformed to simulator coordinates. They appear
C--      in array values in optVar indexed order.
         optVal = values (iopt)
         optVal = optVarLo (iopt) +
     &    values (iopt) * (optVarHi (iopt) - optVarLo (iopt))

C--      Now we have the desired value of the optimization variable
C--      in simulator units. Next, format the number according to
C--      the format specified in subFormat. Add parentheses to the
C--      raw format specifier for this substitution
         write ( optValStr, '(' // subFormat (i) // ')' ) optVal

      end if


C--   The string to be substituted for the mark is now in optValStr.
C--   It appears that optValStr is nicely formatted:
C--     a field width specified by the format statement
C--     the field left justified in the character* string
C--     the number right-justified in the field
C--   That's just what we want. Get its trimmed length.
      optValStrLen = Len_Trim (optValStr)

C--   The output line is made up of three components:
C--    1. All the characters in the original line up to but
C--       not including the {.
C--    2. The trimmed formatted value just calculated
C--    3. All the characters in the original line from the
C--       character after the } to the trimmed end of the string.
      line =
     & line(1:lcb-1) // optValStr (1:optValStrLen) // line(rcb+1:)

C--   We've only promised to do one substitution mark in this routine,
C--   so we're outta here; return success.
      done = .false.
      SubstMark = .true.
      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      real function UserOptValue (iOptVar, nValue)
      implicit none

C--   Given the optimization variable index (1, 2, 3 .. nOptVars)
C--   and the normalized value of the variable, return the value in
C--   the user's coordinate system.

      integer iOptVar
      real nValue

      include 'boss.inc'


      if (iOptVar .gt. nOptVars)  stop

      UserOptValue = optVarLo (ioptVar) +
     &  nValue * (optVarHi (iOptVar) - optVarLo (iOptVar))

      return

      end




C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine TransmitFile (iWorker, fileName)
      implicit none

C--   Transmit a file from the bd subdirectory of the boss to the
C--   wd directory of a specified worker. All errors are
C--   characterized in the diagnostic file, and are fatal.

C--   Worker index (1, 2, 3...)
      integer iWorker

C--   Name of the file to be transmitted (no path required)
      character*(*) fileName

C--   Internal status variable
      integer istat

      include 'boss.inc'

C--   Current worker name
      character*(maxWorkerNameLen) wName

C--   Trimmed length of the current worker name
      integer wNameLen

C--   Source and destination strings
      character*(bigStringLen)  srcString, destString

C--   Trimmed length of the source and destination strings
      integer srcStringLen, destStringLen

C--   External function type
      integer Remove



C--   Get the name of the destination worker
c      write(*,901) iWorker,fileName
c901   format('entry TransmitFile',i6,a20)
      wName = workerName (iWorker)
      wNameLen = Len_Trim(wName)

C--   Path and name of the local file
C--   Example: /home/oc/myFile.dat

      srcString = pathBD // fileName
      srcStringLen = Len_Trim (srcString)

C--   Name of the worker-data directory in the Worker processor,
C--   in scp-compatible format.
C--   Example: opticool@linuxfarm3:wd

      destString =
     &  account // '@' //
     &  wName(1:wNameLen) // ':' //
     &  PathWD
      destStringLen = Len_Trim (destString)

C--   Here goes the actual remote copy command
C--   Example: scp /home/oc/myFile.dat opticool@linuxfarm3:wd
c      write(*,902) pathBD, fileName , workerName(iWorker),
c     & pathWD, relPathWD
c902   format(5(2x,a15))
      call System (
c     &   'rcp ' //
     &   'scp ' //
     &   srcString (1:srcStringLen) // ' ' //
     &   destString (1:destStringLen),
     &   istat)

C--   If the remote-copy fails, it's a fatal error
      if (istat .ne. 0)  then

         write (diagUnit, 800)
     &    wName (1:wNameLen),
     &    srcString (1:srcStringLen),
     &    destString (1:destStringLen),
     &    istat

 800     format (/'FATAL ERROR from routine TransmitFile - - '/
     &         '  Trying to scp local-source to remote-destination'/
     &         '  Destination worker name:      ', a,
     &         '  Local source path/name:       ', a,
     &         '  Remote destination directory: ', a,
     &         '  Status returned from System call: ', i5)

         stop

      end if

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine DeleteFile (fName)
      implicit none

C--   Delete a specified file from the Boss's bd subdirectory.
C--   All errors are reported to the diagnostic file and are fatal.

      include 'boss.inc'

      character*(*) fName

C--   Get around a deficiency in the g77 compiler
      character*(bigStringLen) pathFile

      integer istat

      integer Remove



C--   Construct the pathname/filename string
      pathFile = pathBD // fName(1:Len_Trim(fName)) // '\0'

C--   Remove the file
      istat = Remove (pathFile(1:Len_Trim(pathFile)))

C--   If the remove fails, it's a fatal error
      if (istat .ne. 0)  then

         write (diagUnit, 800)
     &    fName(1:Len_Trim(fName)),
     &    istat

 800     format (/'FATAL ERROR from routine DeleteFile - - '/
     &         '  Trying to delete local boss-data file'/
     &         '  Local filename:              ', a,
     &         '  Status returned from Remove: ', i5)

         stop

      end if

      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine TimeOut (msg, nnn)
      implicit none

C--   Output the current time with an explanatory message.

      include 'boss.inc'

      character*(*) msg
      integer nnn
      integer msgLen
      integer*8 Time8
      character*30  CTime, now



      now = CTime (Time8())

      msgLen = Len_Trim (msg)

      if (nnn .ne. -9999) then

         write (*,800) msg(1:msgLen), nnn
         write (repoUnit,800) msg(1:msgLen), nnn
         write (diagUnit,800) msg(1:msgLen), nnn
 800     format (/'Time Check: ', a, i2)

      else

         write (*,801) msg(1:msgLen)
         write (repoUnit,801) msg(1:msgLen)
         write (diagUnit,801) msg(1:msgLen)
 801     format (/'Time Check: ', a)

      end if

      write (*,802) now
      write (repoUnit,802) now
      write (diagUnit,802) now
 802  format ('Time = ', a/)

      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine AssignWorkload (nParticles)
      implicit none

C--   Given the total number of particles to be propagated, and what
C--   we know about the capabilities of the Workers, assign a
C--   particle propagation count to each living worker.

      integer nParticles

      include 'boss.inc'

      integer nLiveWorkers
      real systemSpeed
      integer iWorker


C--   Find out how many living Workers there are, and what their
C--   total computing power is.

C--   Zero the live-worker counter and the total system speed
      nLiveWorkers = 0
      systemSpeed = 0.

C--   For every worker declared by the user . . .
      do  iWorker = 1, nWorkers

C--      If this worker is alive but idle . . .
         if (workerState (iWorker) .eq. 1)  then

C--         Increment the live worker count
            nLiveWorkers = nLiveWorkers + 1

C--         Add the current computing power to that of the system
            systemSpeed = systemSpeed + workerSpeedCur (iWorker)

         end if

      end do

C--   If there are no living workers, we're stymied
      if (nLiveWorkers .le. 0)  STOP

C--   Allocate the total workload according to capability of the workers
C--   For every living worker . . .
      do  iWorker = 1, nWorkers

         if (workerState (iWorker) .eq. 1) then

C--         Number of particles for this worker is its speed-weighted
C--         share of the whole. Note that workerNPart is an integer
C--         and we don't look after rounding, so the actual total
C--         number of particles thrown may differ from nParticles
C--         by a few.
            workerNPart (iWorker) =
     &       (workerSpeedCur (iWorker) / systemSpeed) * nParticles

         end if

      end do

      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine StartSim (iWorker, randOffset, optValue)
      implicit none

C--   Start a simulation of a specified configuration on a specified
C--   Worker.  The routine returns as soon as the simulation is
C--   initiated.
C--     iWorker:    worker index (1,2,3...nWorkers)
C--     randOffset  random number seed offset (0,1,2,3...99)
C--     optValue    array of nOptVars optimization variable values

      include 'boss.inc'

C--   Worker index (1,2,3...nWorkers)
      integer iWorker

C--   Number of particles for this Worker to propagate
      integer nParts

C--   Offset of the random number seed to be used (0,1,2...)
      integer randOffset

C--   Array of normalized optimization variable values
      real optValue (nOptVars)

C--   True if status is ok
      logical ok

C--   Index over input filename
      integer i

C--   Declare function types
      logical SubstFile



C--   If the worker is not idle, fatal error
      if (workerState (iWorker) .ne. 1)  STOP

C--   Set the random number seed
      workerSeedCur (iWorker) = workerSeedBase (iWorker) + randOffset

C--   For all input files to be transferred . . .
      do  i = 1, nInputFiles

C--      If the modified name is null (no substitution necessary)
         if (Len_Trim (modifFilename (i)) .eq. 0) then

C--         Transmit the file to the worker without modification;
C--         Don't delete this file in the boss.
            call TransmitFile (iWorker, inputFilename (i))

C--      If the modified name is not null, substitution is needed
         else

C--         Perform the substitution
            ok = SubstFile (
     &       pathBD // inputFilename (i),
     &       pathBD // modifFilename (i),
     &       optValue,
     &       iWorker )
            if (.not. ok)  STOP

C--         Transmit the modified file to the worker
            call TransmitFile (iWorker, modifFilename (i))

C--         Delete the boss copy of the modified file
            call DeleteFile (modifFilename (i))

         end if

       end do


C--   Send WORK file to the worker's wd directory; this should
C--   start the simulation on this worker.
      call TransmitFile (iWorker, 'WORK')

C--   Change the worker state to "active"
      workerState (iWorker) = 2
      nWorkersIdle = nWorkersIdle - 1
      nWorkersActive = nWorkersActive + 1

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine WaitResults (iWorker)
      implicit none

C--   Wait until some active worker returns results. When one does,
C--   return its worker index to the caller. This routine blocks
C--   until results are returned; don't call it unless there are
C--   workers in the "active" state.

C--   Index of a worker that has returned results
      integer iWorker

      include 'boss.inc'

      character*2 wIndex
      character*10 resultsFileName

      logical exists



C--   For every worker in the list of workers . . .
 10   do  iWorker = 1, nWorkers

C--      If this Worker is active . . .
         if (workerState (iWorker) .eq. 2) then

C--         Worker index in i2.2 format
            write (wIndex,800) iWorker
 800        format (i2.2)

C--         Name of this Worker's RESULTS file
            resultsFileName = 'RESULTS.' // wIndex

C--         See if the merit output file has arrived yet
            inquire (file = pathBD // resultsFileName, exist=exists)

C--         If this worker has delivered the goods . . .
            if (exists) then

C--            Change the worker's state to "done"
               workerState (iWorker) = 3
               nWorkersActive = nWorkersActive - 1
               nWorkersDone = nWorkersDone + 1

C--            Delete the RESULTS signal file
               call DeleteFile (resultsFileName)

C--            Return the worker index to the caller
               return

            end if

         end if

      end do


C--   No active Worker has returned RESULTS; wait awhile and
C--   try again . . .

      call sleep(1)
      go to 10

      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine FetchResults (iWorker)
      implicit none

      integer iWorker

      include 'boss.inc'

C--   General purpose indexes
      integer i, j

C--   Status returned from file management calls
      integer istat

C--   True if specified file exists
      logical exists

C--   Function type definitions
      logical GetWord

C--   Worker index formatted in i2.2 format (01, 02 ... 10, 11 ...)
      character*2 wIndex

C--   Unit number for reading merit files
      integer merUnit
      parameter (merUnit = recycUnitA)

C--   Actual number of merit variables returned by a simulator
      integer nSimMeritVars

C--   Merit variable names from the simulator's merit file
      character*(maxMeritVarNameLen) simMeritName (maxSimMeritVars)

C--   Merit variable value from the simulator's merit file
      real simMeritValue (maxSimMeritVars)

C--   line of text read from the simulator's merit variable file
      character*(maxSimMeritLineLen) line

C--   Merit name being tested
      character*(maxMeritVarNameLen) meritName

C--   i/o status returned from reads
      integer ios

C--   Index over characters in an input line
      integer index

C--   Merit filename
      character*(bigStringLen) merFileName



C--   Zero counter of merit variables read from the file. Note that
C--   this may include many merit variables not used by the
C--   current optimization.

      nSimMeritVars = 0

C--   Generate file extension (01, 02 . . .) from worker number
      write (wIndex, 800) iWorker
 800  format (i2.2)

C--   Name of this Worker's Merit Output file
      merFileName = 'MeritOut.' // wIndex

C--   Open the merit variable file for this worker
      open (
     &  unit=merUnit,
     &  file = pathBD // merFileName,
     &  status='old',
     &  iostat=ios)

C--   If there's any problem opening the merit file, it's death
      if (ios .ne. 0)  stop

C--   Read-one-line loop . . .
      do  while (.true.)

C--      Don't let the table overflow
         if (nSimMeritVars .ge. maxSimMeritVars) stop

C--      Read one line
         read (merUnit, 805, iostat=ios)  line
 805     format (a)

C--      If it's end-of-file, we're ready to process the merit
C--      variables returned by this worker.
         if (ios .eq. -1)  go to 100

C--      If it's an error, stop
         if (ios .ne. 0)  stop

C--      Increment the line counter
         nSimMeritVars = nSimMeritVars + 1

C--      It's a good line; extract the merit variable name
         index = 1
         if (.not. GetWord (
     &    line, index, simMeritName(nSimMeritVars))) stop

C--      Extract the value of the merit variable
         read (line(index:),*, iostat=ios)
     &    simMeritValue (nSimMeritVars)

C--      Fatal error if extraction fails
         if (ios .ne. 0)  stop

C--   One line decoded; next line . . .
      end do

C--   Close and delete the merit file
      close (merUnit)
      call DeleteFile (merFileName)

C--   The merit file has been read and the results stored in
C--   simMeritName and simMeritValue arrays. Now we must match
C--   these up with the merit variables this optimization is
C--   interested in

 100  continue

C--   Stash the merit values of interest to this optimization
C--   in tables indexed by merit variable index and worker index.

C--   For each merit variable declared in the control file . . .
      do  i = 1, nMeritVars

C--      Get the name of this variable
         meritName = meritVarName (i)

C--      For each merit variable returned by the simulator . . .
         do  j = 1, nSimMeritVars

C--         If it's a name match . . .
            if (simMeritName (j) .eq. meritName) then

C--            Stash the merit variable value
               meritVarValue (i, iWorker) = simMeritValue (j)

C--            No need to search further
               go to 105

            end if

         end do

C--      No match found. That's a fatal error.
         stop

C--      Matched this name; try the next . . .
 105     continue

C--   ...loop over merit variables named in the control file
      end do


      return

      end


C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine ComputeWorkerMerit (iWorker)
      implicit none

      integer iWorker

      include 'boss.inc'


C--   Name of a specific merit variable
      character*(maxMeritVarNameLen) meritName

C--   General purpose index variables
      integer i, j, k

C--   Intermediate results in merit computations
      real x, mNorm

C--   Accumulators for nine merit summations A1..A9
      real meritSum (9)

C--   True if the associated merit sum has been used
      logical meritSumUsed (9)

C--   Merit combination code (M=0, A1=1 . . . A9=9)
      integer comCode

C--   Final figure of merit
      real merit



C--   For every merit variable returned . . .
      do  i = 1, nMeritVars

C--      Transform merit variable value onto (-3, 3)
         x = 3. + 6. *
     &    (meritVarValue (i, iWorker) - meritVarPright (i)) /
     &    (meritVarPright (i) - meritVarPleft (i))

C--      Normalized merit (0 to 1)
         mNorm = (atan (x) / 3.14159265) + 0.5

C--      Actual merit contribution from this merit variable
         meritContrib (i, iWorker) = meritVarMmin (i) +
     &    mNorm * (meritVarMmax (i) - meritVarMmin (i))

      end do


C--   Given the merit contributions, compute the total merit.
C--   First compute the nine subtotals. (A1 ... A9)
      do  k = 1, 9
         meritSum (k) = 0.
         meritSumUsed (k) = .false.
      end do

      do  i = 1, nMeritVars
         comCode = meritVarCombin (i)
         if (comCode .ge. 1 .and. comCode .le. 9) then
            meritSum (comCode) =
     &       meritSum (comCode) + meritContrib (i, iWorker)
            meritSumUsed (comCode) = .true.
         end if
      end do

C--   Now multiply together all of the used sums
      merit = 1.
      do  k = 1, 9
         if (meritSumUsed (k)) merit = merit * meritSum (k)
      end do

C--   Finally, multiply in all of the merit components with M
C--   combination flags
      do  i = 1, nMeritVars
         if (meritVarCombin (i) .eq. 0)
     &    merit = merit * meritContrib (i, iWorker)
      end do


C--   Store this in the table of merits (by worker #)
      meritWorker (iWorker) = merit

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine KillWorkers
      implicit none

      include 'boss.inc'

      integer iWorker


C--   The optimization is done; send kill commands to all the worker
C--   processes, irrespective of their current state. Any results
C--   received from processors after this routine has been called
C--   should be ignored.

C--   For each worker . . .
      do  iWorker = 1, nWorkers

C--      Send DONE file to the worker's wd directory
         call TransmitFile (iWorker, 'DONE')

C--      Change the worker's state to "dead"
         workerState (iWorker) = 0

      end do

C--   All processors are dead
      nWorkersIdle = 0

      nWorkersActive = 0
      nWorkersDone = 0

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM
c
c      subroutine DrawPlots
c      implicit none
c
c      include 'boss.inc'
c
c      integer i
c      integer nColumns, nRows
c
cC---- Output the column labels
c      nColumns = 0
c
cC--   Comments for a human-readable postscript file
c      write (plotUnit, 804)
c 804  format (//'% Column Labels (optimization variable names)'/)
c
cC--   For each optimization variable . . .
c      do  i = 1, nOptVars
c
cC--      If summary simulations are to be performed for this variable
c         if (optVarSumConfigs (i) .gt. 0) then
c
cC--         Increment the column count
c            nColumns = nColumns + 1
c
cC--         Write the column label at the top of the page
c            write (plotUnit, 805) nColumns,
c     &       optVarName(i)(1:Len_Trim(optVarName(i)))
c 805        format ( 1x, i2, ' (', a, ') colLabel' )
c
c         end if
c
c      end do
c
cC--   If no summary simulations were called for, exit now
c      if (nColumns .le. 0)  return
c
c
cC---- Output the row labels
c      nRows = 1
c
cC--   Comments for a human-readable postscript file
c      write (plotUnit, 809)
c 809  format (//'% Row labels (merit variable names)'/)
c
c
cC--   The first row is reserved for the total merit
c      write (plotUnit, 810)
c 810  format (/' 1 (Total Merit) rowLabel')
c
cC--   For each merit variable . . .
c      do  i = 1, nMeritVars
c
cC--      Increment the row count
c         nRows = nRows + 1
c
cC--      Write the row label on the left edge of the page
c         write (plotUnit, 811) nRows,
c     &    meritVarName(i)(1:Len_Trim(meritVarName(i)))
c 811     format ( 1x, i2, ' (', a, ') rowLabel' )
c
c      end do
c
c
cC--   Plot the total merit for the first optimization variable
c      write (plotUnit, 819)
c 819  format (//'% Total Merit, first optimization variable')
c
c      write (plotUnit, 820) 1, 1, optimumNorm(1)
c 820  format (1x, i2, 1x, i2, 1x, ' plotBase'/
c     & ' drawAxes'/
c     & f6.2, ' dashedLine')
c
c      write (plotUnit, 822) optVarLo(1), optVarHi(1)
c 822  format ( ' (
c
c
c
C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM

      subroutine PrintSumConfigs (unit)
      implicit none

      integer unit

      include 'boss.inc'

C--   Index over summary configurations
      integer isc

C--   Index over merit variables
      integer imv

C--   Declare function types
      real UserOptValue



      write (unit,800) nSumConfigs
 800  format (//'Total Number of Summary Configurations: ', i4)


C--   For each summary configuration . . .
      do  isc = 1, nSumConfigs

C--      Print the name and value of the free variable and the
C--      total merit for this configuration
         write (unit,805)
     &    isc,
     &    optVarName (sumFreeIndex (isc)),
     &    UserOptValue (sumFreeIndex (isc), sumFreeValue (isc)),
     &    sumFreeValue (isc),
     &    sumMeritTotal (isc)

 805     format (//
     &    '  Summary configuration number  . . . . . . . ', i12/
     &    '  Free optimization variable name . . . . . . ', a/
     &    '  Free optimization variable user value . . . ', f12.4/
     &    '  Free optimization variable normal value . . ', f12.4/
     &    '  Total Merit for this configuration  . . . . ', f12.4)

C--      For each merit variable . . .
         do  imv = 1, nMeritVars

C--         Print the name of the merit variable, the actual
C--         value of the merit variable, and the merit component
C--         for the merit variable
            write (unit,810)
     &       imv,
     &       meritVarName (imv),
     &       sumMeritVar (isc, imv),
     &       sumMeritComp (isc, imv)

 810        format (/'    Merit variable number . . . . . . ', i12/
     &               '    Merit variable name . . . . . . . ', a/
     &               '    Value of the merit variable . . . ', f12.4/
     &               '    Merit component for the variable  ', f12.4)

         end do

      end do

      return
      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM


      subroutine PrintSumConfigs2 (unit)
      implicit none

      integer unit

      include 'boss.inc'

C--   Index over optimization variables
      integer iov

C--   Number of summary configurations for optimization variable
      integer nsc

C--   Index within summary configuration table
      integer isc

C--   Index within merit variable table
      integer imv

C--   Lower and upper index within summary table for this optVar
      integer ilo, ihi

C--   Declare function types
      real UserOptValue




C--   If there are no summary configurations, return now
      if (nSumConfigs .le. 0)  return

C--   Start at the beginning of the summary table
      ilo = 1


C--   For every optimization variable . . .
      do  iov = 1, nOptVars

C--      If summary configurations were called for . . .
         nsc = optVarSumConfigs (iov)
         if (nsc .gt. 0)  then

C--         Last entry in the summary table for this opt. variable
            ihi = ilo + nsc - 1

C--         Write a header for the optimization variable
            write (unit,800) optVarName (iov)
 800        format (//'Summary Configurations for Optimization ',
     &       'Variable ', a)

C--         Write a header for the total-merit section
            write (unit,805)
 805        format (/'  Total Configuration Merit')

C--         Write the merits for these summary configurations
            do  isc = ilo, ihi

               write (unit,810)
     &          UserOptValue (sumFreeIndex (isc), sumFreeValue (isc)),
     &          sumFreeValue (isc),
     &          sumMeritTotal (isc)

 810           format ('    ', f12.4, ' ( ', f7.4, ' ): ', f12.4)

            end do

C--         For each merit variable . . .
            do  imv = 1, nMeritVars

C--            Write a header for the merit variable
               write (unit,815) meritVarName (imv)
 815           format (/'  Merit Variable Value and Merit Component ',
     &          'for Merit Variable ', a)

C--            For each summary configuration . . .
               do  isc = ilo, ihi

                  write (unit,820)
     &             UserOptValue (sumFreeIndex(isc), sumFreeValue(isc)),
     &             sumFreeValue (isc),
     &             sumMeritVar (isc, imv),
     &             sumMeritComp (isc, imv)

 820              format ('    ', f12.4, ' ( ', f7.4, ' ): ', 2f12.4)

               end do

            end do

C--         Step ahead in the summary table
            ilo = ihi + 1

         end if

      end do

      return

      end

C MWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWMWM
