Recipe: Simulated Annealing v1.01

created by LociOiling

Profile

• Name: Simulated Annealing v1.01
• ID: 48981
• Shared with: Public
• Parent: Simulated Annealing v1.0
• Children: None
• Created on: April 28, 2014 at 06:12 AM UTC
• Updated on: April 28, 2014 at 06:12 AM UTC
• Description:

  Performs a great many operations on a protein while using Simulated Annealing algorithm (see wikipedia) to decide which ones to accept and which to reject. Reference to FuseBurstIterateCi function removed as temporary fix for missing code.

Code

--------------------------------------------------------------------
--           Code for Simulated Annealing algorithm           --
--------------------------------------------------------------------
-- Performs "modifications" of a protein and accepts or rejects them depending on "temperature"
-- and on how much they improved the score (if + improvement, yes; if not too bad and temperature
-- is still high, then probably yes). Temperature starts very high, drops pretty quickly.
--
--  On a Win7 box with ~50 AAs, a 100 step run takes 30 min to 1 hour; with ~100 AAs, 90 min is closer
--  Plan runs accordingly.
--
-- ALLCAP variables are intended for users to modify as desired:
-- 1. In early stages (or if a protein is "stuck" hard), one might want lots of restarts with not 
--     many steps before a restart 
--    (eg STEPSPERRUN = 20, NUMRESTARTS = 50, REHEAT_NOT_RESTART = false)
-- 2. Later stages might prefer to increase STEPSPERRUN, decrease NUMRESTARTS, and
--     set REHEAT_NOT_RESTART to true.
-- 3. Mutations: if a puzzle allows mutations, but a segment should not be mutated, it can
--     be added to the NOMUTATELIST. If a puzzle allows mutations, but one doesn't want
--     them to happen at all, set FORBIDMUTATION to true.
--
-- TODO: 
-- 1.  Include a band-two-sheets-together banding function and a band-two-cysteines together one, too.
-- 2.  Go over: Quake Rebuilders, idealize (not micro), deep rebuilders, alternative compressors - in particular,
--      idealize and rebuild are both pretty primitive. Provide some combo-actions involving them
-- 3.  Create "WorkCores" - something that does lots of core-compression runs and then short SAs
--      on any that are interesting. Save them in Quicksave slots so SA or later scripts can find them
-- 4.  Upgrade the gui so users can directly choose actions (eg 1 bander/freezer + rebuild + 1 wiggle; WorkCores)
-- FARTHER FUTURE
-- 5.  Collect stats about how each modifier and each wiggler performed (time and score). Report them and try
--      modifying probabilities of choosing those routines on the fly. (also modify relative probs of small/medium/
--      large-scale actions)
-- 6. Study the Fuse-wigglers, sort them out, decide which to use/keep. Play with freeze-wigglers, make sure 
--     Wiggle and LocalWiggle are mixing well (both All and Selected for each, with selected either a range or
--     an "interesting" collection - such as one that expands and shrinks its selected list based on what moved), 
-- 7. Does CleanupWiggle really need to use iterative wiggles? Simpler ones would be quicker and might
--     do well enough for scoring needs. Perhaps instead iterative wiggles should happen *between* SA runs, or
--     only if SA decides to accept a neighbor?
--
-- Code ideas and/or thefts come from: 
--   (tlaloc) tlaloc Random Tug 4.00
--   (spvincent) Helix Twister 1.0, Loop rebuild 5.0 (NC), Local Quake 1.0
--   (susume2) Cut and Wiggle Everything v0.1
--   (MurloW) Shock v0.1, Fracture v1.6
--   (MurloW, drjr) Idealize by 3+
--   (porkythepundit) NC Kichen KaboodlEnz 14, Porky's Cement Crusher 1.0 (ideas for fuses)
--   (rav3n_pl) Rav3n_pl GAB BiS v2.0.1 w30 [NC], Rav3n_pl Voids Killer v0.5 NC, Rav3n_pl Push v3, LS Quake v1.1
--   (StackOverflow) Freeze Burst v1.0
--   (tlaloc, rav3n_pl, Seagat2011, thom001) ST - Glycine Hinge REPOST
--   (drjr) drjr - RandomizeR 1.20, Idealize by 3
--
--  revision history
--  v1.01   2014/04/27 LociOiling - commented reference to missing FuseBurstIterateCi
--------------------------------------------------------------------------------

SA_STEPSPERRUN = 100.0
SA_NUMSTARTS = 4
SA_MINGAIN = 1.0
MF_STEPSPERRUN  = 50.0
MF_FAILSBEFORERESTART = 30
MARCHWIGGLE_MINGAIN = 2.0
TARG_STEPSPERRUN = 10
TARG_FAILSBEFORERESTART = 10

CORE_PULLINGCI = 0.70

-- values for controlling how mutations are handled (under-tested)
FORBIDMUTATION = false        -- this is for performing local "spot" changes
MUTATE_NOT_SHAKE = false  -- this is for choosing between shake and mutate
NOMUTATELIST = {0}
PROB_MUTATESIMILAR = 0.75

-- values for controlling wiggle choices and timing
QUICKWIGGLE_ITERS = 2
SLOWWIGGLE_ITERS = 2
FINALWIGGLE_ITERS = 8

-- values for controlling behavior of slow wiggle that's done after an acceptable neighbor is found
PROB_SLOWWIGGLE_ALL = 0.60
PROB_SLOWWIGGLE_LOCAL = 0.30
-- PROB_SLOWWIGGLE_CUTFUSE = 0.10  -- 1.00 - PROB_SLOWWIGGLE_ALL - PROB_SLOWWIGGLE_LOCAL

QUICK_IT_WIGGLE_GAIN  = 0.50    -- for fuses that want iterative wiggles even in "quickwiggle" time
SLOW_IT_WIGGLE_GAIN    = 0.10   -- for cleanup wiggle actions at each step
FINALRUN_WIGGLE_GAIN = 0.001  -- for final "drag out points" activities

-- values for choosing the "scale" of modification to try
PROB_TARGETED = 0.05
PROB_LOCAL = 0.50
PROB_MIDSIZE = 0.40
-- "PROB_LARGESCALE" is 1.00 - PROB_LOCAL - PROB_MIDSIZE - PROB_TARGETED

-- values used in idealize
USE_CUTS = true
FORBID_UNHEALED_CUTS = true

-- some oddballs that might want setting
USENORMALSCORE  = true  -- exploration puzzles would set false
DEBUGRUN                = true   -- mostly goes with DebugPrint, but could get more use later

-------------- The following are not generally useful to change --------------
SA_FAILSBEFORERESTART   = 100  -- doesn't seem useful so far
REHEAT_NOT_RESTART = true  -- mid/late-game use strongly prefers true

-- values for adjusting rebuild parameters, probably not wise for callers to adjust
INNER_TUBE_RADIUS  = 12.0  -- used by rebuilds; 9.0 has been seen in some scripts
OUTER_TUBE_RADIUS = 18.0  -- used by rebuilds


-------------------------------------------------------------------------------------------------------------
------------ VALUES USERS SHOULD NOT BE TOUCHING ---------------------------------
-------------------------------------------------------------------------------------------------------------
-- QUICKSAVE slots used herein
QS_Start = 1
QS_Run  = 2
QS_Best = 3

QS_SATemp = 91
QS_NeighborTemp = 92
QS_MFTemp = 93    -- might merge with SATemp as a "high-level" controller function storage place
QS_ShortUseTemp = 98   -- for use by low-level subroutines that just want to stash state for short time
QS_ShortUseTemp2 = 99 -- for use by low-level subroutines that just want to stash state for short time

PuzzleAllowsMutate = false
PuzzleHasLockedSegs = false
InitialScore = 0.0
InitialSelectedIdxes = {}
InitialSecondaryIdxes = {}
InitialFrozenTable = {}
InitialBandCount = 0
InitialEnabledBandTable = {} -- not using because disabled band between segs prevents new band there
IsDesignerPuzzle = false
EndCalled = false
IdealSimpleHelixDistance4 = 6.3  -- for skips of 4
IdealSimpleHelixDistance5 = 8.6  -- for skips of 5
StartTime = os.time()
helixStarts = {}
helixEnds = {}
sheetStarts = {}
sheetEnds = {}
loopStarts = {}
loopEnds = {}
NonLockedSegList = {}
InitialClashImportance = behavior.GetClashImportance()
segCt = structure.GetCount()
whichAlg = ""


function DebugPrint( str )
    if DEBUGRUN then print( str ) end
end

----------------------------------------------------------------
-- MATHY STUFF
----------------------------------------------------------------
function TrimNum(val)
    return val - val % 0.001 
end

function seedRandom()
    seed=os.time()/math.abs(getScore())   -- Do not copy blindly: get getScore() or replace with current.GetEnergyScore()
    seed=seed%0.001
    seed=1/seed
    while seed<10000000 do seed=seed*1000 end
    seed=seed-seed%1
    -- DebugPrint("Seed is: "..seed)
    math.randomseed(seed)
end

function random( n1,n2, forceFloat ) --random function returns int or float depends on input vars
    if forceFloat == nil then forceFloat = false end
    if n1==nil then 
        return math.random()
    else
        if n2==nil then
            if n1%1==0 then -- can't test for "forceFloat", so caller must beware
                return math.random(n1) --integer
            else
                return math.random() * n1 --float
            end
        else
            if n1%1==0 and n2%1==0 and not forceFloat then
                return math.random( n1, n2 ) --integer between 
            else
                return math.random() * (n2 - n1) + n1 --float between
            end
        end
    end
end

function randomDice( ctDice, minValPerDie, maxValPerDie )  -- distro that prefers "middle" values
    total = 0
    for i=1, ctDice do
        total = total + random( minValPerDie, maxValPerDie )
    end
    return total
end

function randomThetaPhi()
    return   math.acos( random( -1.0, 1.0, true ) ),  random( 2 * math.pi )
end

function randomizeIndexList( idxList )
    for i=1, #idxList do
        j = random( #idxList )
        if j ~= i then
            idxList[i], idxList[j] = idxList[j], idxList[i]
        end
    end
end

function randomSeg()
    return random( segCt )
end

function randomMovableSeg()
    if not PuzzleHasLockedSegs then 
        return randomSeg() 
    end

    return NonLockedSegList[ random( #NonLockedSegList ) ]
end

function randomLowScoringSeg( ctIndicesToChooseAmong, subscore )
    local idxList = {}
    getWorstScoringAAs( idxList, ctIndicesToChooseAmong, subscore )
    return idxList[ random(  #idxList ) ]
end

function Coprime( n )
   -- find the highest prime < 70% of "n" that is coprime with "n"
    local primes = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101}
    for i = #primes, 1, -1 do
        if primes[ i ] < 0.70*n and n % primes[ i ] ~= 0 then
            return primes[ i ]
        end
    end
    return 1
end

function GetEuclideanDistanceVector( scmat1, scmat2, vecOut )
    score = 0.0

    for j=1, #scmat1[1] do
        vecOut[j] = 0.0
    end

    for i=1, segCt do
        for j=1, #scmat1[i] do
            vecOut[j] = vecOut[j] + (scmat2[i][j] - scmat1[i][j]) * (scmat2[i][j] - scmat1[i][j])
        end
    end

    for j=1, #vecOut do
        vecOut[j] = math.sqrt( vecOut[j] )
    end
end

----------------------------------------------------------------
-- INFORMATION GATHERING
----------------------------------------------------------------
function getScore( wantRB )
    if wantRB == nil then wantRB = false end
    local s=0.0
    if not USENORMALSCORE then
        if wantRB then s = recentbest.GetEnergyScore()
        else s=current.GetEnergyScore()
        end
    else
        if wantRB then s = recentbest.GetScore()
        else s=current.GetScore()
        end
    end
    return s
end

function ConstructScoreMatrix( scmat )
    -- just what *are* the various subscores? Various sources name quite a pile of possibilities
    -- leaving out "reference" (unmodifiable anyway), "loctotal", "total", "ligand"
    for i=1, segCt do
        scmat[i] = {}
        scmat[i][1]   =  current.GetSegmentEnergySubscore( i, "backbone" )
        scmat[i][2]   =  current.GetSegmentEnergySubscore( i, "sidechain" )
        scmat[i][3]   =  current.GetSegmentEnergySubscore( i, "clashing" )
        scmat[i][4]   =  current.GetSegmentEnergySubscore( i, "packing" ) 
        scmat[i][5]   =  current.GetSegmentEnergySubscore( i, "density" )
        scmat[i][6]   =  current.GetSegmentEnergySubscore( i, "hiding" )
        scmat[i][7]   =  current.GetSegmentEnergySubscore( i, "bonding" )
        scmat[i][8]   =  current.GetSegmentEnergySubscore( i, "disulfides" )
        scmat[i][9]   =  current.GetSegmentEnergySubscore( i, "ideality" )
        scmat[i][10]   =  current.GetSegmentEnergySubscore( i, "pairwise" )
        scmat[i][11] =  current.GetSegmentEnergySubscore( i, "other" )
    end
end

function checkAla()
    local ala=0
    for n = 1, segCt do
        segType = structure.GetAminoAcid(n)
        if segType == 'a' or segType == 'g' then
            ala=ala+1
        end
    end
    return ala == segCt
end

function IsPuzzleMutable()
   for i=1,segCt do
      if structure.IsMutable(i) then return true end
   end
   return false
end
function isNoMutate( idx )
   if FORBIDMUTATION then return true end
   if not PuzzleAllowsMutate then return true end
   if not structure.IsMutable(idx) then return true end

   for i=1,#NOMUTATELIST do
       if idx == NOMUTATELIST[i] then return true end
   end
   return false
end

function FindAllMovableSegs()
    for i=1,segCt do
        if structure.IsLocked(i) then 
            PuzzleHasLockedSegs = true
        else
            NonLockedSegList[ #NonLockedSegList + 1] = i
        end
    end
    return false
end

function isMovableSeg( seg1 )
   if seg1 == nil or seg1 < 1 or seg1 > segCt then return false end

   local BB, SC = freeze.IsFrozen( seg1 )
   local sl = structure.IsLocked( seg1 )
   return not (BB or SC or sl)
end

function isAllMovableRange(seg1, seg2)
    if seg1 == nil or seg2 == nil then return false end
    if seg1 > seg2 then seg1, seg2 = seg2, seg1 end
    if seg1 < 1 or seg2 > segCt then return false end

    for i=seg1, seg2 do
        if not isMovableSeg( i ) then return false end
    end
    return true
end

function RemoveCutAtIndex( idx )
    if idx == nil or idx < 1 or idx > segCt then return true end
    local sc = getScore()
    structure.DeleteCut( idx )
    return getScore() ~= sc  -- if no score change at all, then cut didn't happen!
end

------------- REGION-IDENTIFICATION FUNCTIONS ------------
function getHelices()
    inside = false
    for i=1, segCt do
        if ( structure.GetSecondaryStructure(i) == "H" ) then
            if ( inside == false ) then
                inside = true
                helixStarts[ #helixStarts + 1] = i
            end
        elseif ( inside == true ) then
            inside = false
            helixEnds[ #helixEnds + 1] = i - 1
        end -- if ( "H" ) elseif ( within )
    end -- for (segCt)

    -- deal with "last seg is helix"
    if ( inside == true ) then
        helixEnds[ #helixEnds + 1] = segCt
    end
end

function getSheets()
    inside = false
    for i=1, segCt do
        if ( structure.GetSecondaryStructure(i) == "E" ) then
            if ( inside == false ) then
                inside = true
                sheetStarts[ #sheetStarts + 1 ] = i
            end
        elseif ( inside == true ) then
            inside = false
            sheetEnds[ #sheetEnds + 1 ] = i - 1
        end -- if/else 'E'
    end -- for (segCt)

    -- deal with "last seg is sheet"
    if ( inside == true ) then
        sheetEnds[ #sheetEnds + 1 ] = segCt
    end
end

function getLoops()
    inside = false
    for i=1, segCt do
        if ( structure.GetSecondaryStructure(i) == "L" ) then
            if ( inside == false ) then
                inside = true
                loopStarts[ #loopStarts + 1 ] = i
            end
        elseif ( inside == true ) then
            inside = false
            loopEnds[ #loopEnds + 1 ] = i - 1
        end -- if/else 'L'
    end -- for (segCt)

    -- deal with "last seg is loop"
    if ( inside == true ) then
        loopEnds[ #loopEnds + 1 ] = segCt
    end
end

function getRegionStartAndEnd( idx ) 
    s = idx
    e = idx
    local typeIdx = structure.GetSecondaryStructure( idx )
    if idx > 1 then
        s = idx - 1
        while s >= 1 do
            if structure.GetSecondaryStructure( s ) ~= typeIdx then break end
            s = s - 1
        end
        s = s + 1
    end
     if idx < segCt then
        e = idx + 1
        while e <= segCt do
            if structure.GetSecondaryStructure( e ) ~= typeIdx then break end
            e = e + 1
        end
        e = e - 1
    end
    return s, e
end

-- uses getRegionStartAndEnd, but tries to tune "loop" regions a bit
function getRegionForOperation( idx )
    -- try to find an acceptable region (really bizarre choices of SS can give trouble)
    sIdx, eIdx = getRegionStartAndEnd( idx )

    -- Check: if our entire region is only 1 or 2 segs long, then see about expanding it
    local regionType = structure.GetSecondaryStructure( idx )
    if eIdx - sIdx <= 1 then 
        -- small range. Try to expand it...
        local sIdx2, eIdx2
        if sIdx > 2 then
            if structure.GetSecondaryStructure(sIdx - 1) == regionType then
                sIdx2, eIdx2 = getRegionStartAndEnd( sIdx - 1)
                sIdx = sIdx2
            end
         end
        if eIdx < segCt - 1 then
            if structure.GetSecondaryStructure( eIdx + 1 ) == regionType then
                sIdx2, eIdx2 = getRegionStartAndEnd( eIdx + 1 )
                eIdx = eIdx2
            end
        end -- END try to fix small range
        if eIdx - sIdx <= 1 then return idx,idx end  -- failed
    end

    -- Check: if loop and more than 30-long, then only take part of it
    if regionType == 'L' and eIdx - sIdx > 30 then
        -- trim it; it is too long
        if idx - sIdx < 15 then
            eIdx = idx + 15  -- trim the other, longer side
        elseif eIdx - idx < 15 then
            sIdx = idx - 15  -- trim the other, longer side
        else
            -- more than 15 segs slop on both sides!
            -- try to be random and "medium-long": 3d5 on each side should do
            sIdx = idx - randomDice(3, 1, 5)
            eIdx = idx + randomDice(3, 1, 5)
        end
    end

    return sIdx, eIdx
end


---------- SORTING FOR DISTANCES AND SCORES ------------
function sortAllByScore(segs, scores)
    table.sort(segs, function(n1,n2) return scores[n1] < scores[n2] end)
end

------------------- SCORING FUNCTIONS ------------------
-- returns all aas in order of score (if subscore is nil then does GetSegmentEnergyScore)
function getWorstScoringAAs( idxList, maxWanted, subscore )
    local scoreList = {}
    local idxes = {}
    for i=1, segCt do
       idxes[i] = i
    end
    for i=1, segCt do
        if subscore == nil then
            scoreList[i] = current.GetSegmentEnergyScore(i)
        else
            scoreList[i] = current.GetSegmentEnergySubscore( i, subscore )
        end
    end
    sortAllByScore( idxes, scoreList )
    for i = 1, math.min( maxWanted, segCt ) do
        idxList[i] = idxes[i]
    end
end
function getBestScoringAAs( idxList, maxWanted, subscore )
    local scoreList = {}
    local idxes = {}
    for i=1, segCt do
       idxes[i] = i
    end
    for i=1, segCt do
        if subscore == nil then
            scoreList[i] = current.GetSegmentEnergyScore(i)
        else
            scoreList[i] = current.GetSegmentEnergySubscore( i, subscore )
        end
    end
    sortAllByScore( idxes, scoreList )
    for i = 1, math.min( maxWanted, segCt ) do
        idxList[i] = idxes[segCt - i + 1]
    end
end
function isBadScorer(idx, mustBeAtLeastThisBad)
    local idxList = {}
    getWorstScoringAAs( idxList, mustBeAtLeastThisBad)
    for i=1, #idxList do
        if idxList[i] == idx then return true end
    end
    return false
end

------------------- DISTANCE FUNCTIONS ------------------
function computeDistanceSums(distList)
    for i=1, segCt do
        distList[i] = 0.0 
    end
    for i=1, segCt-1 do
        for j=i+1, segCt do
            distList[i] = distList[i] + structure.GetDistance( i, j ) / segCt
            distList[j] = distList[j] + structure.GetDistance( i, j ) / segCt
        end
    end
end
function getCentralAAs( idxList )
    local distList = {}
    for i=1, segCt do
       idxList[i] = i
    end
    computeDistanceSums( distList )
    sortAllByScore( idxList, distList )
end

-- returns all hydrophobics in order of "centralness"
function getCentralHydrophobics(idxList, ctWanted)
    local idxListInternal = {}
    getCentralAAs(idxListInternal)
    local ct = 1
    for i=1, #idxListInternal do
        if structure.IsHydrophobic( idxListInternal[ i ] ) then
            idxList[ ct ] = idxListInternal[ i ]
            if ct == ctWanted then return end
            ct = ct + 1
        end
    end
end
function isCentralHydrophobic(idx, mustBeAtLeastThisGood)
    local idxList = {}
    getCentralHydrophobics( idxList, mustBeAtLeastThisGood)
    for i=1, #idxList do
        if idxList[i] == idx then return true end
    end
    return false
end

function getSegsWithinSphere( idxList, centerIdx, radius, includeCenter )
    if includeCenter then
        idxList[ #idxList + 1 ] = centerIdx
    end
    for i=1, segCt do
        if i ~= centerIdx and structure.GetDistance( centerIdx, i ) <= radius then
            idxList[ #idxList + 1 ] = i
        end
    end
end

function getSegsOutsideSphere( idxList, centerIdx, radius )
    for i=1, segCt do
        if i ~= centerIdx and structure.GetDistance( centerIdx, i ) > radius then
            idxList[ #idxList + 1 ] = i
        end
    end
end

function getSegsWithinRangeOfTube( idxList, startSeg, endSeg, radius )
    -- if a seg is "near" to one of the ones in our [startSeg, endSeg] range, then we want it
    for i=1, segCt do
        for j=startSeg, endSeg do
            if structure.GetDistance( i, j ) < radius then
                idxList[ #idxList + 1] = i
                break
            end
        end  -- for j
    end      -- for i
end

function IsLocalMinimum( segCenter,  segCheck)
    if segCheck == segCenter then return false end  -- exclude selfies
    ckDist = structure.GetDistance(segCenter, segCheck)
    if segCheck > 1 then 
        ckDist1 = structure.GetDistance( segCenter, segCheck - 1 )
        if ckDist1 < ckDist then return false end
        if segCheck > 2 then 
            ckDist2 = structure.GetDistance( segCenter, segCheck - 2 )
            if ckDist2 < ckDist then return false end
        end
    end
    if segCheck < segCt then
        ckDist1 = structure.GetDistance( segCenter, segCheck + 1 )
        if ckDist1 < ckDist then return false end
        if segCheck < segCt - 1 then 
            ckDist2 = structure.GetDistance( segCenter, segCheck + 2 )
            if ckDist2 < ckDist then return false end
        end
    end
    return true
end

function IsLocalMaximum( segCenter,  segCheck)
    if segCheck == segCenter then return false end  -- exclude selfies
    ckDist = structure.GetDistance(segCenter, segCheck)
    if segCheck > 1 then 
        ckDist1 = structure.GetDistance( segCenter, segCheck - 1 )
        if ckDist1 > ckDist then return false end
        if segCheck > 2 then 
            ckDist2 = structure.GetDistance( segCenter, segCheck - 2 )
            if ckDist2 > ckDist then return false end
        end
    end
    if segCheck < segCt then
        ckDist1 = structure.GetDistance( segCenter, segCheck + 1 )
        if ckDist1 > ckDist then return false end
        if segCheck < segCt - 1 then 
            ckDist2 = structure.GetDistance( segCenter, segCheck + 2 )
            if ckDist2 > ckDist then return false end
        end
    end
    return true
end

------------------------------------------------------------------------------
--
---- WIGGLER FUNCTIONS
--
------------------------------------------------------------------------------

-- simple shake/mutate function
function ShakeOrMutate( time, idxList)
    if time == nil then time = 1 end
    if idxList == nil then
        if MUTATE_NOT_SHAKE and not FORBIDMUTATION then
            structure.MutateSidechainsAll( time )
        else
            if not IsDesignerPuzzle then structure.ShakeSidechainsAll( time ) end
        end
    else
        selection.DeselectAll( )
        for i=1, #idxList do selection.Select( i ) end
        if MUTATE_NOT_SHAKE and not FORBIDMUTATION then
            structure.MutateSidechainsSelected ( time )
        else
            if not IsDesignerPuzzle then structure.ShakeSidechainsSelected ( time ) end
        end
        selection.DeselectAll( )
    end
end

------------------------------------------------------------------------------
-- QUICK WIGGLERS
------------------------------------------------------------------------------
function WiggleRange( time, startSeg, endSeg, addFreeze, pureLocal)
    local maxSeg = structure.GetCount()
    local lastSeg = math.min( endSeg, maxSeg )
    local firstSeg = math.max( 1, startSeg )
    if lastSeg < firstSeg then return end -- nothing to do!

    -- set up for the local wiggle
    selection.DeselectAll()
    if addFreeze then
        if firstSeg > 1 then freeze.Freeze( firstSeg - 1, true, true ) end
        if lastSeg < maxSeg then freeze.Freeze( lastSeg + 1, true, true ) end
    end

    -- actually perform the local wiggle
    selection.SelectRange( firstSeg, lastSeg )
    if pureLocal then
        structure.LocalWiggleSelected( time, true, true )
    else
        structure.WiggleSelected( time, true, true )
    end

    -- clean up what we did
    if addFreeze then
        if firstSeg > 1 then freeze.Unfreeze( firstSeg - 1, true, true ) end
        if lastSeg < maxSeg then freeze.Unfreeze( lastSeg + 1, true, true ) end
    end
end

function WiggleLocalByChunk(time, startSeg, endSeg, chunkSize, passCount, addFreeze)
    currOffset = 0
    for i=1, passCount do
        if (passCount == 1 or chunkSize == 1) then currOffset = 0
        else currOffset = random(chunkSize - 1)
        end

        ctBlocks = (endSeg - startSeg) / chunkSize
        for idx = 0,ctBlocks do
            startIdx = startSeg + currOffset + chunkSize*idx
            WiggleRange(time, startIdx, startIdx + chunkSize, addFreeze, true)
        end
    end
end

function WiggleWalk(time, doForward, doPureLocal)
    for i=1,segCt do
        if doForward then idx = i
        else idx = segCt + 1 - i
        end
        WiggleRange(time, idx, idx, false, doPureLocal)
    end
end

------------------------------------------------------------------------------
-- ITERATIVE SIMPLE WIGGLERS
------------------------------------------------------------------------------
function IterativeWiggleAll( time, minGain, doBackbone, doSidechain )
    local lastScore = getScore( )
    local gain = minGain
    while gain >= minGain do
        structure.WiggleAll( time, doBackbone, doSidechain )  -- do's ok for nil
        gain = getScore( ) - lastScore
        lastScore = getScore( )
    end
end

function IterativeWiggleList( idxList, time, minGain, doBackbone, doSidechain )
    selection.DeselectAll( )
    for i=1, #idxList do selection.Select( i ) end

    local lastScore = getScore( )
    local gain = minGain
    while gain >= minGain do
        structure.WiggleSelected( time, doBackbone, doSidechain )  -- do's ok for nil
        gain = getScore( ) - lastScore
        lastScore = getScore( )
    end
    selection.DeselectAll( )
end

function IterativeWiggleLocalByChunk(time, startSeg, endSeg, minGain, chunkSize, addFreeze)
    currOffset = 0
    local lastScore = getScore( )
    gain = minGain
    while gain >= minGain do
        if chunkSize == 1 then currOffset = 0
        else currOffset = random(chunkSize - 1)
        end

        ctBlocks = (endSeg - startSeg) / chunkSize
        for idx = 0,ctBlocks do
            startIdx = startSeg + currOffset + chunkSize*idx
            WiggleRange(time, startIdx, startIdx + chunkSize, addFreeze, true)
             if addFreeze then freeze.UnfreezeAll() end
        end
        local currScore = getScore( )
        gain = currScore - lastScore
        lastScore = currScore
    end
end

-- give it a list of segments to perform its actions upon
function IterativeWiggleShakeWiggle( idxList, time )
    IterativeWiggleList( idxList, time, QUICK_IT_WIGGLE_GAIN )
    ShakeOrMutate( 1, idxList )
    IterativeWiggleList( idxList, time, QUICK_IT_WIGGLE_GAIN )
end

function qStabWiggle( time, doExtraWork )
    print("QStabWiggle")
    if doExtraWork == nil then doExtraWork = false end
    local ciInit = behavior.GetClashImportance()

    behavior.SetClashImportance( 0.10)
    ShakeOrMutate( 1 )

    behavior.SetClashImportance( 1.0)
    if doExtraWork then
        IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
        ShakeOrMutate( 1 )
    end
    IterativeWiggleAll( time, 10*QUICK_IT_WIGGLE_GAIN, true, true ) 

    behavior.SetClashImportance( ciInit )
end

--------------------------------------------------
-- FUSE WIGGLERS
--------------------------------------------------
function WiggleShakeSimpleFuse( ci, time )
    behavior.SetClashImportance( ci )
    structure.WiggleAll( time )  -- instead of IterativeWiggleAll...
    ShakeOrMutate( time/2 )
end
function WiggleWiggleSimpleFuse( ci1, ci2, shorttime, longtime )
    behavior.SetClashImportance( ci1 )
    structure.WiggleAll( shorttime )
    behavior.SetClashImportance( ci2 )
    structure.WiggleAll( longtime )  -- "long" wiggle
end
function ShakeWiggleFuse( ci1, ci2, time )  -- aka Fuze1
    behavior.SetClashImportance( ci1 )
    ShakeOrMutate( time/2 )
    behavior.SetClashImportance( ci2 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
end
function WiggleWiggleWiggleFuse( ci1, ci2, time ) -- aka Fuze2
    behavior.SetClashImportance( ci1 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
    behavior.SetClashImportance( 1.00 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
    behavior.SetClashImportance( ci2 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
end
function WiggleShakeWiggleFuse( time )  -- aka FuzeEnd
    behavior.SetClashImportance( 1.00 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
    ShakeOrMutate( 1 )
    IterativeWiggleAll( time, QUICK_IT_WIGGLE_GAIN, true, true )
end

function QuickFuseBurst ( basetime )
    local temp = behavior.GetClashImportance( )
    WiggleShakeSimpleFuse( 0.01, basetime )
    WiggleShakeSimpleFuse( 0.50, basetime )
    WiggleShakeSimpleFuse( 1.00, basetime*2 )
    behavior.SetClashImportance ( temp )
end

function BlueFuse( time )
    local temp = behavior.GetClashImportance( )
    -- another blue has 0.05-shake(1), 1-wiggleall(8), 
    -- 0.07-shake(1), 1-wiggleall(8)+check, 
    -- 0.30-wiggleall(1), 1-wiggleall(8)+check
    ShakeWiggleFuse( 0.02, 1.00, time )
    ShakeWiggleFuse( 0.25, 1.00, time )
    ShakeWiggleFuse( 0.55, 1.00, time )
    ShakeWiggleFuse( 1.00, 1.00, time )
    behavior.SetClashImportance ( temp )
end

function PinkFuse( time )
    local temp = behavior.GetClashImportance( )
    ShakeWiggleFuse( 0.10, 0.70, time ) 
    WiggleShakeWiggleFuse( time )
    ShakeWiggleFuse( 0.30, 0.60, time )
    WiggleShakeWiggleFuse( time )
    WiggleWiggleWiggleFuse( 0.50, 0.70, time )
    WiggleShakeWiggleFuse( time )
    WiggleWiggleWiggleFuse( 0.70, 0.50, time )
    WiggleShakeWiggleFuse( time )
    behavior.SetClashImportance ( temp )
end

function CementCrusherFuse( time )
    local temp = behavior.GetClashImportance( )
    CIlow = {0.98, 0.89, 0.55, 0.34, 0.21, 0.13, 0.08, 0.05, 0.03, 0.02, 0.01}
    for i=1, #CIlow do
        WiggleWiggleSimpleFuse( CIlow[i], 1.00, 1, time )
    end
    behavior.SetClashImportance ( temp )
end

------------------------------------------------------------------------------
-- SLOW WIGGLERS
------------------------------------------------------------------------------
-- big pile of Fuse attempts; if keepGoing is true, can be very slow
function FuseWiggleNicely( time, keepGoing )
    if keepGoing == nil then keepGoing = false end
    local ciInit = behavior.GetClashImportance()
    recentbest.Save()
    rbScore = getScore( true )
    print("FuseWiggleNicely with CI="..ciInit)

    -- gentlemen, start your fuzing!

    ShakeWiggleFuse( 0.3, 0.6, time )
    WiggleShakeWiggleFuse( time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    WiggleWiggleWiggleFuse( 0.5, 0.7, time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    ShakeWiggleFuse( 0.05, 1, time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    WiggleWiggleWiggleFuse( 0.7, 0.5, time )
    WiggleShakeWiggleFuse( time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    ShakeWiggleFuse( 0.07, 1, time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    WiggleWiggleWiggleFuse( 0.3, 1, time ) 
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    -- getting desperate for points: try BlueFuse
    BlueFuse( time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    -- one more try: PinkFuse
    PinkFuse( time )
    if getScore( true ) < rbScore then 
        recentbest.Restore() 
    elseif not keepGoing then 
        behavior.SetClashImportance( ciInit )
        return 
    end

    -- and a truly final last try: CementCrusher
    CementCrusherFuse( time )
    if getScore( true ) < rbScore then recentbest.Restore() end

    behavior.SetClashImportance( ciInit )
end

function PerformFreezeBurst( baseTime )
    local frozenCt = structure.GetCount () - 2
    while ( frozenCt > 1 ) do
        recentbest.Save()
        local idxList = {}
        getBestScoringAAs( idxList, frozenCt )
        for i = 1, frozenCt do
            freeze.Freeze ( idxList[i] , true , true )
        end
        QuickFuseBurst ( baseTime )
        freeze.UnfreezeAll ()
        recentbest.Restore()
        frozenCt = frozenCt / 2
    end
    freeze.UnfreezeAll ()
end

function PerformCutAndWiggleAction( basetime, shift, offset )
    if shift == nil then shift = 3 end
    if offset == nil then offset = 1 else offset = 1 + (offset % shift) end
    local temp = behavior.GetClashImportance( )

    for i = offset, segCt-1, shift do
        structure.InsertCut(i)
    end
    WiggleWiggleSimpleFuse( 0.30, 1.0, basetime, 2*basetime )

    for i = 1, segCt, shift do -- take the opportunity to kill any "lying around" cuts
        RemoveCutAtIndex( i )
    end
    WiggleWiggleSimpleFuse( 0.30, 1.0, basetime, 4*basetime )

    behavior.SetClashImportance ( temp )
end

-- good, but expensive
function IterativeCutAndWiggle( basetime, minGain, shift )
    local offset = random( shift )
    local score = getScore()
    local gain = minGain
    repeat
        recentbest.Save()
        PerformCutAndWiggleAction( basetime, shift, offset )
        local newscore = getScore()
        if newscore < score then 
            recentbest.Restore()
            gain = 0.0
        else 
            gain = newscore - score
            score = newscore
        end
        offset = 1 + ( offset % shift )  -- offset + 1 works, but ashamed to use it
    until gain < minGain
end

-- very expensive
function FullIterativeCutAndWiggle( basetime, minGain, maxShift, randomizeShifts )
    repeat
        local gain = 0.0
        local score = getScore()
        local newscore = score
        local shift = maxShift
        if randomizeShifts then shift = random(3,maxShift) end
        for offset = 1, shift do
            recentbest.Save()
            PerformCutAndWiggleAction( basetime, shift, offset )
            newscore = getScore()
            if newscore < score then 
                recentbest.Restore()
            end
        end  -- FOR offset
        gain = newscore - score
    until gain < minGain
end

--------------------------------------------------
--BAND-CHANGING WIGGLERS
--------------------------------------------------
-- returns the actual amount of change that occurred (absolute value)
function BandedWiggleUntilEnoughChange( time, scoreThreshold, pullingCI )
    if pullingCI == nil then pullingCI = CORE_PULLINGCI end

    temp = behavior.GetClashImportance()
    behavior.SetClashImportance( pullingCI )
    deltaPoints = 0.0
    ss=getScore()

    for str=0.30, 1.00, 0.11 do --search enough band strength to move
        for i = InitialBandCount+1, band.GetCount() do
            band.SetStrength(i, str)
        end
        structure.WiggleAll( time, true, false)
        if math.abs( ss-getScore() ) > scoreThreshold then 
            -- we've changed the structure "enough", send it off 
            break 
        end
    end
    behavior.SetClashImportance( temp )
    return math.abs( ss - getScore() )
end

-- returns the actual amount of change that occurred (absolute value)
function BandedWiggleUntilTargetMet( time, scoreDeltaTarget, maxTries, pullingCI )
    save.Quicksave( QS_ShortUseTemp )
    save.Quicksave( QS_ShortUseTemp2 )
    local temp = behavior.GetClashImportance()
    if pullingCI == nil then pullingCI = CORE_PULLINGCI end
    behavior.SetClashImportance( pullingCI )

    local s1=getScore()
    local str=0.50
    local tries = 0
    local s2best = 99999
    while tries < maxTries do
        save.Quickload( QS_ShortUseTemp)
        for i = InitialBandCount + 1, band.GetCount() do band.SetStrength(i, str) end
        structure.WiggleAll( time )

        local s2 = math.abs( getScore() - s1 )
        if s2 > scoreDeltaTarget * 1.2 then 
            str = str * ( 0.85 + random( 0.10 ) )
        elseif s2 < scoreDeltaTarget * 0.8 then 
            str = str + ( 1.0 - str ) * ( 0.05 + random( 0.1 ) )
        else 
            return math.abs( getScore() - s1 )   -- we're good enough
        end
        if math.abs( s2 - scoreDeltaTarget ) < s2best then
            save.Quicksave( QS_ShortUseTemp2 )
            s2best = math.abs( s2 - scoreDeltaTarget )
        end

        tries = tries + 1
    end
    save.Quickload( QS_ShortUseTemp2 )  -- the best we could come up with...
    behavior.SetClashImportance( temp )
    return math.abs( getScore() - s1 )  -- returns the actual delta in case caller cares
end

function SimpleBandedWiggleWrapper( time, minChangeWanted, pullingCI )
    local scoreThreshold = math.min( minChangeWanted, getScore() / 500.0)
    BandedWiggleUntilEnoughChange( time, scoreThreshold, pullingCI )
    DelBands()
    qStabWiggle( time, random() < 0.50 )  -- a bit of wiggle afterwards to help settle down
end

function SimpleWiggleBeforeAndAfter( time, startIdx, endIdx, probQstab )
    -- make a change. Perhaps measure the change for "enoughness"?
    local wantGlobalWiggle = random() < 0.50
    if random() < probQstab then
        qStabWiggle( time, random() < 0.50 ) -- 2nd param is a "work harder" one
    else
        SimpleLocalWiggle( time, startIdx, endIdx, wantGlobalWiggle, random(4) )
    end

    freeze.UnfreezeAll()  -- if we froze anything, we should clean that up now
    DelBands()  -- if we put in bands, they should go now

    if random() < probQstab then
        qStabWiggle( time, random() < 0.50 ) -- 2nd param is a "work harder" one
    else
        SimpleLocalWiggle( time, startIdx, endIdx, not wantGlobalWiggle, random(4,8) ) -- range bigger than previous
    end
end

-------------------------------------------
-- NORMAL WIGGLE WRAPUP --
-------------------------------------------
function LocalWiggleNicely( time, skipCt, passCount )
    print( "WiggleLocalByChunk with skip="..skipCt )
    WiggleLocalByChunk( time, 2, segCt, skipCt, passCount, false)
    ShakeOrMutate( 1 )
    structure.WiggleAll(2,true,true)
end

function WiggleWalkNicely( time )
    local r = random(4)
    print( "WiggleWalk with r="..r )
    WiggleWalk( time, r < 3, r == 1 or r == 3 )
    ShakeOrMutate( 1 )
    structure.WiggleAll( time ,true, true )
end

function WiggleZlbNicely( time, strength, skipLower, skipUpper )
    print( "WiggleZlbNicely with strength="..strength )
    if skipLower == nil then skipLower = segCt+1 end
    if skipUpper == nil then skipUpper = 0 end
    ZLB(strength, skipLower, skipUpper)
    ShakeOrMutate( 2 )
    structure.WiggleAll(time, true, true)
end

function SimpleLocalWiggle( time, idxLo, idxHi, pickGlobalWiggle, wiggleExpansion )
    -- caller gives it a core range (ok for Lo == Hi). Routine expands range "randomly".

    if wiggleExpansion == nil then
        wiggleExpansion = 1 + random( math.floor( segCt / 6 ))
    end
    local startIdx = math.max( 1, idxLo - wiggleExpansion )
    local endIdx = math.min( segCt, idxHi + wiggleExpansion )
    print( "SimpleLocalWiggle from ".. startIdx .." to ".. endIdx )
    WiggleRange( time, startIdx, endIdx, false, pickGlobalWiggle )
end

--------------------------------------------------------------
--
---- MUTATOR FUNCTIONS
-- TODO: Consider adding a "very similar" mutator
--------------------------------------------------------------
aaCount = 20
aaPhilCount = 9
aaPhobCount = 11
function AA(idx)
    -- hydrophilics
    if idx == 1 then return "r"
    elseif idx == 2 then return "s"
    elseif idx == 3 then return "t"
    elseif idx == 4 then return "n"
    elseif idx == 5 then return "d"
    elseif idx == 6 then return "q"
    elseif idx == 7 then return "e"
    elseif idx == 8 then return "h"
    elseif idx == 9 then return "k"
    -- hydrophobics
    elseif idx == 10 then return "g"
    elseif idx == 11 then return "a"
    elseif idx == 12 then return "c"
    elseif idx == 13 then return "v"
    elseif idx == 14 then return "l"
    elseif idx == 15 then return "i"
    elseif idx == 16 then return "m"
    elseif idx == 17 then return "p"
    elseif idx == 18 then return "f"
    elseif idx == 19 then return "y"
    else return "w" end      -- idx == 20
end

function AAphilic(idx)
    -- hydrophilics
    if idx == 1 then return "r"
    elseif idx == 2 then return "s"
    elseif idx == 3 then return "t"
    elseif idx == 4 then return "n"
    elseif idx == 5 then return "d"
    elseif idx == 6 then return "q"
    elseif idx == 7 then return "e"
    elseif idx == 8 then return "h"
    else return "k"      -- idx == 9
    end
end

function AAphobic(idx)
    -- hydrophobics
    if idx == 1 then return "g"
    elseif idx == 2 then return "a"
    elseif idx == 3 then return "c"
    elseif idx == 4 then return "v"
    elseif idx == 5 then return "l"
    elseif idx == 6 then return "i"
    elseif idx == 7 then return "m"
    elseif idx == 8 then return "p"
    elseif idx == 9 then return "f"
    elseif idx == 10 then return "y"
    else return "w"      -- idx == 11
    end 
end

-- this one doesn't care what amino used to be at idx
function MutateToAny( idx )
    local aa = structure.GetAminoAcid( idx )
    local aanew = aa
    local i = 0
    while ( aanew == aa ) do
        i = random( aaCount )
        aanew = AA( i )
    end
    print( "mutating " .. aa .. " at " .. idx .. " to " .. aanew )
    structure.SetAminoAcid( idx, aanew )
    return true
end

-- this one preserves hydrophob/hydrophil
function MutateToSimilar( idx )
    local aa = structure.GetAminoAcid( idx )
    local i = 0
    local aanew = aa
    while (aanew == aa) do
        if structure.IsHydrophobic( idx ) then
            i = random( aaPhobCount )
            aanew = AAphobic( i )
        else
            i = random( aaPhilCount )
            aanew = AAphilic( i )
        end
    end
    print( "mutating " .. aa .. " at " .. idx .. " to " .. aanew )
    structure.SetAminoAcid( idx, aanew )
    return true
end


----------------------------------------------------------------------
--
---- BANDER FUNCTIONS
--
----------------------------------------------------------------------
function DelBands()
    --for i=InitialBandCount+1, band.GetCount() do
    --    band.Delete( i )
    --end
    band.DeleteAll()
end

function BandBetweenSegsWithParameters( seg1, seg2, strength, goalLength )
    -- caller can set goalLength without setting strength by providing strength = 0.0
    if not ( isMovableSeg( seg1) or isMovableSeg( seg2) ) then return false end

    local bIdx = band.AddBetweenSegments( seg1, seg2 )
    if bIdx ~= band.GetCount( ) then 
        DebugPrint( "failed to add band from "..seg1.." to "..seg2)
        return false 
    end

    if goalLength ~= nil then 
        band.SetGoalLength( bIdx, goalLength ) 
    end
    if strength ~= nil and strength > 0.0 then 
        band.SetStrength( bIdx, strength )
    end
    return true
end

function BandInSpaceWithParameters( seg, segFini, segInit, rho, theta, phi, strength, goalLength )
    -- caller can set goalLength without setting strength by providing strength = 0.0
    if not ( isMovableSeg( seg) ) then return false end

    local bIdx = band.Add( seg, segFini, segInit, rho, theta, phi )
    if bIdx ~= band.GetCount( ) then return false end

    if goalLength ~= nil then 
        band.SetGoalLength( bIdx, goalLength ) 
    end
    if strength ~= nil and strength ~= 0.0 then 
        band.SetStrength( bIdx, strength ) 
    end
    return true
end

----------------------------------------------------------------------
-- SMALL-SCALE BANDERS
-- consider putting a "desired length" (or desired fraction of dist between segs)
-- try making a 2-nearest-cystein joiner
----------------------------------------------------------------------
function PutBandInSpace( idx, maxRho, strength )
    local idx2, idx3
    if idx < segCt then
        idx2 = idx + 1
    else
        idx2 = idx - 2
    end
    if idx > 1 then
        idx3 = idx - 1
    else
        idx3 = idx + 2
    end

    -- random point in sphere of radius maxRho
    local theta, phi = randomThetaPhi( )
    local rho = (maxRho * random()^(1/3)) + 0.001
    return BandInSpaceWithParameters( idx, idx2, idx3, rho, theta, phi, strength )
end

function PutBandToRandomSeg( idx, strength )
    needNew = true
    local failedTries = 0
    while ( needNew and failedTries < 20 ) do
        idx2 = random( segCt )  -- ok if this one isn't movable (we assume idx is movable)
        if idx2 > idx + 1 or idx2 < idx - 1 then 
            needNew = false 
            break
        else
            failedTries = failedTries + 1
        end
    end
    if not needNew then
        return BandBetweenSegsWithParameters( idx, idx2, strength )
    end
    return false
end

function PutGlycineHingeBands( glycineIdx, strength )
    if glycineIdx == 1 or glycineIdx == segCt then return false end
    changeSucceeded = false
    selection.DeselectAll ()

   -- construct the hinge (freeze before and after the glycine)
    local sR, eR = getRegionStartAndEnd( glycineIdx + 1 )
    if eR < segCt and eR - glycineIdx < 3 then eR = eR + 1 end
    selection.SelectRange ( glycineIdx + 1 , eR )

    local sL, eL = getRegionStartAndEnd( glycineIdx - 1 )
    if eL > 1 and glycineIdx - eL < 3 then eL = eL - 1 end
    selection.SelectRange ( sL, glycineIdx - 1 )
    freeze.FreezeSelected ( true, true )
    selection.DeselectAll ()

    -- put bands connecting the two rods 
    local maxLen = math.min( eR - sR + 1, eL - sL + 1 )
    if maxLen > 4 then maxLen = 4
    elseif maxLen < 1 then return changeFailed  -- should only happen if glycine at an end
    end
    for i=1, maxLen do
        if BandBetweenSegsWithParameters( math.max(1, glycineIdx - i), math.min(segCt, glycineIdx + i), strength ) then 
            changeSucceeded = true 
        end
    end
    return changeSucceeded
end

----------------------------------------------------------------------
-- MID-SCALE BANDERS
-- still want a 2-sheet-joiner collection of bands
----------------------------------------------------------------------

function PutPushPullBands( badSeg, bandsWanted, center, wantPushNotPull )
    changeSucceeded = false

    local okSegList = {}
    local ct = 1
    local SC = structure.GetDistance( badSeg, center )
    local SR

    for i=1, segCt do
        local RC=structure.GetDistance( i, center )
        SR=structure.GetDistance( i, badSeg )
        if SR<15
            and RC<SC-3
            and math.abs( i - badSeg ) > 5
            and math.abs( i - center ) > 3
        then
            okSegList[ct] = i
            ct = ct + 1
        end
    end
    randomizeIndexList( okSegList )

    local bandsAdded = 0
    for i=1, #okSegList do
        SR=structure.GetDistance( i, badSeg )
        if wantPushNotPull then
            SR=SR + 4
        else
            SR=SR - 4
        end
        if SR>20 then SR=20 end
        if SR<3 then SR=3 end

        if BandBetweenSegsWithParameters( badSeg, okSegList[i], 0, SR ) then 
            changeSucceeded = true 
            bandsAdded = bandsAdded + 1
        end
        if bandsAdded >= bandsWanted then break end
    end --  FOR(i)
    return changeSucceeded
end

function PutVoidCrusherBands( idx ) --make band if found void in area of segment "idx"
    changeSucceeded = false
    local t={} -- store possible segments

    for b=1,segCt do --test all segments
        if math.abs(idx-b) >= 15 then    -- else idx,b too near each other count-wise
            local ab = structure.GetDistance(idx, b)
            if ab > 10.0 then --no void if less
                local void=true

                for c=1,segCt do --searching for any segment between them
                    if c ~= idx and c ~= b then   -- don't want c = idx or b
                        local ac = structure.GetDistance(idx,c)
                        local bc = structure.GetDistance(b,c)
                        if ac<ab and bc<ab and ac>4 and bc>4 then
                            if ac+bc<ab+1.5
                                then void=false break --no void there for sure
                            end
                        end
                    end
                end  -- END for c

                if void==true then 
                     t[#t+1]={idx,b}
                end
            end  -- END if idx and b are spatially at least 10 apart
        end -- END if idx and b are count-wise at least 15 apart
    end -- LOOP over candidate b's

    if #t>0 then
        for i=1,#t do
            local s=t[i][1]
            local e=t[i][2]
            local d=structure.GetDistance( s, e )
            d = d - 3.0  -- try to compress this much...
            if d<3 then d=3 end
            if d>20 then d=20 end
            if BandBetweenSegsWithParameters( s, e, 0, d ) then changeSucceeded = true end
        end
    end

    return changeSucceeded
end

-- mostly intended for helices 
function PutUniformDistanceBands( startSeg, endSeg, shift, strength )
    changeSucceeded = false
    sumDist = 0
    for i=startSeg, endSeg-shift do
        sumDist = sumDist + structure.GetDistance( i, i+shift )
    end
    length = sumDist/( segCt - shift )

    for i=startSeg,endSeg-shift do
        if BandBetweenSegsWithParameters( i, i+shift, strength, length ) then changeSucceeded = true end
    end
     return changeSucceeded
end

-- only intended for helices 
-- this might change to allow other offsets such as 3, 6, 7. (larger isn't forbidden)
-- Any future call for double helices would need different ideal distances
function PutHelixFixerBands( startSeg, endSeg, strength, do4, do5 )
    changeSucceeded = false

    if do4 then
        for i=startSeg,endSeg-4 do
            if BandBetweenSegsWithParameters( i, i+4, strength, IdealSimpleHelixDistance4 ) then changeSucceeded = true end
        end
    end
    if do5 then
        for i=startSeg,endSeg-5 do
            if BandBetweenSegsWithParameters( i, i+5, strength, IdealSimpleHelixDistance5 ) then changeSucceeded = true end
        end
    end
    return changeSucceeded
end

-- mostly intended for helices
function PutCompressorBands( startSeg, endSeg, shift, strength, stretch )
    changeSucceeded = false
    for i=startSeg,endSeg-shift do
        dist = structure.GetDistance( i, i+shift )
        if BandBetweenSegsWithParameters( i, i+shift, strength, dist*stretch ) then changeSucceeded = true end
     end
     return changeSucceeded
end

function PutRotatorBands( startSeg, endSeg, length, strength, wantClockwise )
    changeSucceeded = false
    for i = math.max( startSeg , 2 ), math.min( endSeg , segCt - 1 ) do
        if wantClockwise then
            if BandInSpaceWithParameters( i, i + 1, i - 1, length, math.pi/2, math.rad ( 315 ), strength ) then
                changeSucceeded = true
            end
        else
            if BandInSpaceWithParameters( i, i - 1, i + 1, length, math.pi/2, math.rad ( 315 ), strength ) then
                changeSucceeded = true
            end
        end
    end
    return changeSucceeded
end

function PutPusherBands( startSeg, endSeg, length, strength, doForward )
    changeSucceeded = false
    if startSeg > 1 then 
        init = startSeg - 1
    else
        init = startSeg
    end
    if endSeg < segCt then
        fini = endSeg + 1
    else
        fini = endSeg
    end
    if not doForward then init, fini = fini, init end
    for i = init + 1, fini - 1 do
        if BandInSpaceWithParameters( i, fini, init, length, math.pi/2, 0.0, strength ) then
            changeSucceeded = true
        end
    end
    return changeSucceeded
end

function PutMoverBands( startSeg, endSeg, length, strength )
    changeSucceeded = false
    if startSeg > 1 then 
        init = startSeg - 1
    else
        init = startSeg
    end
    if endSeg < segCt then
       fini = endSeg + 1
    else
       fini = endSeg
    end

    theta, phi = randomThetaPhi( )
    for i = init + 1, fini - 1 do
        if BandInSpaceWithParameters( i, fini, init, length, theta, phi, strength ) then
            changeSucceeded = true
        end
    end
    return changeSucceeded
end

-- TODO consider putting more bands in "between" segments, with lengths that increase
-- as they head towards the midpoint
function PutFlexerBands( startSeg, endSeg, length, strength )
    if endSeg - startSeg < 5 then return false end
    changeSucceeded = false
    if startSeg > 1 then
        init = startSeg - 1
        first = startSeg
    else
        init = startSeg
        first = startSeg + 1
    end
    if endSeg < segCt then
       fini = endSeg + 1
       last = endSeg
    else
       fini = endSeg
       last = endSeg - 1
    end
    midpt = (startSeg + endSeg + 1) / 2

    -- initial band points inwards
    if BandInSpaceWithParameters( first, fini, init, 1.0, math.pi/2.0, 0.0, strength ) then
        changeSucceeded = true
    end

    --final band points inwards, too
    if BandInSpaceWithParameters( last, init, fini, 1.0, math.pi/2.0, 0.0, strength ) then
        changeSucceeded = true
    end

    -- middle one is longer and points perpendicularly to init-fini line. Consider
    -- replacing with ramping-up lengths over all "interior" segments.
    theta, phi = randomThetaPhi( )    -- ignoring phi, because we just want to use 0
    if BandInSpaceWithParameters( midpt, fini, init, length, theta, 0.0, strength ) then
        changeSucceeded = true
    end

    return changeSucceeded
end

function PutMultipleZLBandsInSpace( maxBandCt, maxStrength )
    changeSucceeded = false
    local bandCt = 0
    local ctBandsToMake = random( maxBandCt )
    while bandCt < ctBandsToMake do
        local seg = randomMovableSeg()
        local strength = random( 0, maxStrength, true ) + 0.001
        local maxRho = 0.30
        if PutBandInSpace( seg, maxRho, strength ) then
            changeSucceeded = true
            bandCt = bandCt + 1
        end
    end
    return changeSucceeded
end

function PutMultipleBandsBetween( maxBandCt, maxStrength )
    local ctBandsToMake = random( maxBandCt )

    local bandCt = 0
    local ctTries = 0
    while bandCt < ctBandsToMake and ctTries < 1000 do
        local idx = randomMovableSeg()
        local strength = random( 0.0, maxStrength, true )
        if PutBandToRandomSeg( idx, strength ) then
            bandCt = bandCt + 1
            ctTries = ctTries + 1
        end
    end
    return (bandCt == ctBandsToMake)
end

function PutMultipleBandsBothTypes( maxBandBet, maxBIS, maxBBstrength, maxBISstrength )
    changeSucceeded = false
    local ctBetween = random( maxBandBet )
    local ctBis = random( maxBIS )

    if PutMultipleBandsBetween( ctBetween, maxBBstrength ) then
        changeSucceeded = true
    end
    if PutMultipleZLBandsInSpace( ctBis, maxBISstrength ) then
        changeSucceeded = true
    end

    return changeSucceeded
end

-- bandToNearest false means "band to farthest"
function PutManyBandsToSeg( idx, baseStrength, maxBandCount, wantStrengthScaled, bandToNearest, wantFreeze )
    changeSucceeded = false
    local segsToBand = {}
    local idxList = {}
    local distList = {}
    local totalDistance = 0.0
    local dist = 0.0

    -- we will construct segsToBand first, summing totalDistance as we go
    for i = 3, segCt - 3 do
        if i ~= idx and isMovableSeg( i ) then
            if bandToNearest then
                if IsLocalMinimum( idx, i ) then
                    segsToBand[#segsToBand + 1] = i
                    totalDistance = totalDistance + structure.GetDistance( idx, i )
                end
            else
                if IsLocalMaximum( idx, i ) then
                    segsToBand[#segsToBand + 1] = i
                    totalDistance = totalDistance + structure.GetDistance( idx, i )
                end
            end
        end
    end

    -- now sort them so we can take only the "best" segments...
    for i=1, #segsToBand do idxList[i] = i end
    for i=1, #segsToBand do distList[i] = structure.GetDistance( idx, i ) end

    sortAllByScore( idxList, distList )

    -- now put the bands, preferring the "nearest" ones
    for i=1, math.min( #segsToBand, maxBandCount ) do
        local dist = structure.GetDistance( idx, segsToBand[ idxList[ i ] ] )
        if wantStrengthScaled then
            strength = baseStrength * (totalDistance / #segsToBand) / dist  -- WHY??? wouldn't goallength be better?
        else
            strength =  baseStrength
        end
        if BandBetweenSegsWithParameters( idx, segsToBand[ idxList[ i ] ], strength ) then
            changeSucceeded = true
        end
    end

    if changeSucceeded and wantFreeze then 
        local segStart, segEnd = getRegionForOperation( idx )
        selection.DeselectAll( )
        selection.SelectRange( segStart, segEnd )
        freeze.FreezeSelected( true, true )
        selection.DeselectAll( )
    end

    return changeSucceeded
end

----------------------------------------------------------------------
-- LARGE-SCALE BANDERS
----------------------------------------------------------------------
function ZLB(strength, skipLower, skipUpper)
    if skipLower == nil then skipLower = segCt + 1 end
    if skipUpper == nil then skipUpper = 0 end

    for i=2, segCt-1 do 
        if i < skipLower or i > skipUpper then
            BandInSpaceWithParameters( i, i+1, i-1, 0.01, 0.0, 0.0, strength, 0.01)
        end
    end
    if skipLower > 1 then
        BandInSpaceWithParameters( 1, 2, 3, 0.01, 0.0, 0.0, strength, 0.01)
    end
    if skipUpper < segCt then
        BandInSpaceWithParameters( segCt, segCt-1, segCt-2, 0.01, 0.0, 0.0, strength, 0.01)
    end
end

-- for compressing towards the center (or the opposite)
function PutCoreCompressBands(coreList, skipCt, strength, compressFactor)
    local maxToBand = #coreList
    changeSucceeded = false

    -- first put bands between all of our central segs
    -- notice: there's no protection from some of them being adjacent, but it's ok
    -- if one or more of these bands fail to be created.
    for i=1, maxToBand-1 do
        for j=i+1, maxToBand do
            if coreList[i] < coreList[j] - 1 or coreList[i] > coreList[j] + 1 then
                length = compressFactor * structure.GetDistance( coreList[ i ], coreList[ j ] )
                if BandBetweenSegsWithParameters( coreList[ i ], coreList[ j ], strength, length ) then
                    changeSucceeded = true
                end
            end
        end
    end

    -- now put bands between the "rest" of the protein and our chosen ones
    for i=1, segCt, skipCt do
        for j=1, maxToBand do
            length = compressFactor * structure.GetDistance(  i, coreList[ j ] )
            if BandBetweenSegsWithParameters(  i, coreList[ j ], strength, length ) then
                changeSucceeded = true
            end
        end
    end
    return changeSucceeded
end

function PutBandsEveryXToEveryY( Xrange, Xinit, Yrange, Yinit, strength )
    changeSucceeded = false
    for i=Xinit, segCt, Xrange do
        for j=Yinit, segCt, Yrange do
            if (i < j - 1 or i > j + 1) then
                if BandBetweenSegsWithParameters(  i, j, strength ) then
                    changeSucceeded = true
                end
            end
        end
    end
    return changeSucceeded
end

----------------------------------------------------------------------
--
---- IDEALIZERS AND REBUILDERS
--
----------------------------------------------------------------------
----------------------------------------------------------------------
-- IDEALIZERS
----------------------------------------------------------------------
-- cut points are problematic. For now, I simply keep fighting
-- them off every chance I get.
function IdealizeInSmallRegion( idx, time, useCuts )
    if useCuts == nil then useCuts = true end
    local idxList = {}
    changeFailed = false

    if useCuts then
        if idx - 2 >= 1 then 
            print( "  putting cut at ".. idx - 2 )
            structure.InsertCut( idx - 2 ) 
        end
        if idx + 1 <= segCt then 
            print( "  putting cut at ".. idx + 1 )
            structure.InsertCut( idx + 1 ) 
        end
    end

    selection.DeselectAll()
    getSegsWithinSphere( idxList, idx, 4.0, true )
    for i=1, #idxList do selection.Select( idxList[i] ) end
    structure.IdealizeSelected()
    selection.DeselectAll()

    if useCuts then
        if idx - 3 >= 1 then
            BandBetweenSegsWithParameters(  idx-3, idx, 1.0, structure.GetDistance( idx-3, idx ) )
        end

        if idx+3 <= segCt then
            BandBetweenSegsWithParameters(  idx, idx + 3, 1.0, structure.GetDistance( idx, idx+3 ) )
        end

        -- get some wiggling happening
        WiggleZlbNicely( time, 2.0, math.max( 1, idx-2 ), math.min( segCt, idx+2 ) )

        -- now put the protein back in its proper shape
        if idx - 2 >= 1 then
            if not RemoveCutAtIndex( idx - 2 ) and FORBID_UNHEALED_CUTS then changeFailed = true end
        end
        if changeFailed == false and idx + 1 <= segCt then 
            if not RemoveCutAtIndex( idx + 1 ) and FORBID_UNHEALED_CUTS then changeFailed = true end
        end
    else
        -- since we didn't put in any cuts, we can try simple wiggling to see how we did
        changeFailed = false
        WiggleRange( time, idx - 2, idx + 1, false, random() < 0.50 )
    end

    DelBands()
    return not changeFailed
end

----------------------------------------------------------------------
-- REBUILDERS
----------------------------------------------------------------------
-- experiments: freeze before and after start-end, cut before and after start-end
function RebuildInSmallRegion( startSeg, endSeg, maxTries, time )
    changeSucceeded = false
    local oldScoreMat = {}
    local score = getScore( )
    ConstructScoreMatrix( oldScoreMat )

    selection.DeselectAll()
    selection.SelectRange ( startSeg, endSeg )
    save.Quicksave( QS_ShortUseTemp )
    save.Quicksave( QS_ShortUseTemp2 )

    bestRescore = -99999.0
    for i=1, maxTries do
        save.Quickload( QS_ShortUseTemp)
        structure.RebuildSelected( i ) -- not only up to maxTries times, but up to maxTries iterations per time
        if math.abs( getScore( ) - score ) >= 0.01 then
            ShakeOrMutate( time/2, idxListInner )
            structure.WiggleSelected( time )
            rescore = getScore()
            if rescore > score then
                save.Quicksave( QS_ShortUseTemp )
                save.Quicksave( QS_ShortUseTemp2 )
                score = rescore
                bestRescore = rescore
                changeSucceeded = true    -- we did succeed in changing something
            else
                if IsNeighborDifferentFromOld( oldScoreMat, 0.1 ) and rescore > bestRescore then
                    save.Quicksave( QS_ShortUseTemp2 )
                    bestRescore = rescore
                    if score - rescore < 100.0 then changeSucceeded = true end
                end
            end
        end
    end
    selection.DeselectAll( )
    save.Quickload( QS_ShortUseTemp2 )
    return changeSucceeded
end

----------------------------------------------------------------------
--
---- OPERATION-CHOOSING FUNCTIONS
--
----------------------------------------------------------------------

function PickRebuildOperation( startSeg, endSeg, maxRebuildCount )
    -- TODO: values of radii of inner/outer spheres default to 12.0 and 18.0. consider experimenting
    print( "Rebuild small region for "..startSeg.. " to ".. endSeg )
    local time = QUICKWIGGLE_ITERS
    changeSucceeded = false
    local idxListInner = {}
    local idxListOuter = {}
    getSegsWithinRangeOfTube( idxListInner, startSeg, endSeg, INNER_TUBE_RADIUS )
    getSegsWithinRangeOfTube( idxListOuter, startSeg, endSeg, OUTER_TUBE_RADIUS )

    changeSucceeded = RebuildInSmallRegion( startSeg, endSeg, maxRebuildCount, time )

    if changeSucceeded then
        IterativeWiggleShakeWiggle( idxListOuter, time  )
    end
    return changeSucceeded
end

function PickPointRebuildOperation( idx, maxRebuildCount )
    local startIdx = math.max( 1, idx - random(3) )
    local endIdx = math.min( segCt, idx + random(3) )
    return PickRebuildOperation( startIdx, endIdx, maxRebuildCount )
end

function PickRegionBander( startSeg, endSeg )
    local bStren = 1.0
    local bLen = 5.0
    local time = QUICKWIGGLE_ITERS
    local allowNonSheetOperation = structure.GetSecondaryStructure( startSeg ) ~= 'E'
    local allowHelixOnlyOperation = structure.GetSecondaryStructure( startSeg ) == 'H'
    local regionLength = endSeg - startSeg + 1

    local ctActions = 3
    if regionLength > 5 then
        -- only allow Flexer and the helix ones if there at least 6 segs to affect
        ctActions = ctActions + 1   -- Flexer is acceptable
        if allowNonSheetOperation then ctActions = ctActions + 2 end
        if allowHelixOnlyOperation then ctActions = ctActions + 1 end
    end

    local r = random(ctActions)
    if r == 1 then
        print( "Rotate-bands from " .. startSeg.. " to "..endSeg )
        changeSucceeded = PutRotatorBands( startSeg, endSeg, bLen, bStren, random() < 0.50 )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    elseif r == 2 then
        print( "PushForOrBack-bands from " .. startSeg.. " to "..endSeg)
        changeSucceeded = PutPusherBands( startSeg, endSeg, bLen, bStren, random() < 0.50 )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    elseif r == 3 then
        print( "PushDirection-bands from " .. startSeg.. " to "..endSeg)
        changeSucceeded = PutMoverBands( startSeg, endSeg, bLen, bStren )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    elseif r == 4 then
        print( "Flex-bands from " .. startSeg.. " to "..endSeg)
        changeSucceeded = PutFlexerBands( startSeg, endSeg, bLen, bStren )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    elseif r == 5 then 
        local hShift = 4
        if random() < 0.25 then 
            hShift = math.min(regionLength - 2, randomDice(4, 1, 2)) -- a number between 4 and 8, preferring 6ish
        end
        print( "UDist-bands from " .. startSeg.. " to "..endSeg)
        changeSucceeded = PutUniformDistanceBands( startSeg, endSeg, hShift, bStren )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    elseif r == 6 then
        local hShift = 4
        if random() < 0.25 then 
            hShift = math.min(regionLength - 2, randomDice(4, 1, 2)) -- a number between 4 and 8, preferring 6ish
        end
        local hComp = random( 0.75, 1.25 )
        print( "C/Ex-bands from " .. startSeg.. " to "..endSeg)
        changeSucceeded = PutCompressorBands( startSeg, endSeg, hShift, bStren, hComp )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    else --if r == 7 then
        print( "HelixFixer-bands from " .. startSeg.. " to "..endSeg)
        do4 = random() < 0.5
        do5 = random() < 0.5
        changeSucceeded = PutHelixFixerBands( startSeg, endSeg, bStren, do4, do5 )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, startSeg, endSeg, 0.20 ) end
    end

    DelBands()
    return changeSucceeded
end

function PickCentralBander( maxCenters, maxStrength )
    local time = QUICKWIGGLE_ITERS
    local coreSegs = {}
    getCentralHydrophobics( coreSegs, maxCenters)

    local compression = random( 0.75, 1.25, true )
    local strength = random( 0.0, maxStrength, true )
    local skipCt = math.floor( segCt / random(1, 6) )

    print("Central-bands with "..#coreSegs.." centers and expansion="..compression)
    changeSucceeded = PutCoreCompressBands( coreSegs, skipCt, strength, compression )
    if changeSucceeded then SimpleBandedWiggleWrapper( time, 100.0, CORE_PULLINGCI ) end

    DelBands()
    return changeSucceeded
end

function PickEveryXToEveryYBanding( maxStrength )
    local time = QUICKWIGGLE_ITERS
    local Xinterval = random( 3, math.floor( segCt / 4 ) )
    local Yinterval = random( 3, math.floor( segCt / 4 ))
    local Xinit = random( math.floor( segCt / 8 ) )
    local Yinit = random( math.floor( segCt / 8 ) )
    local strength = random( 0.0, maxStrength, true )

    print("XY bands every "..Xinterval.." to every "..Yinterval)
    changeSucceeded = PutBandsEveryXToEveryY( Xinterval , Xinit, Yinterval, Yinit, strength )
    if changeSucceeded then SimpleBandedWiggleWrapper( time, 100.0, CORE_PULLINGCI ) end

    DelBands()
    return changeSucceeded
end

function PickPushPullAction( badSeg, wantPushNotPull )
    local segList = {}
    local time = QUICKWIGGLE_ITERS
    changeSucceeded = false
    DelBands()

    getSegsOutsideSphere( segList, badSeg, 8 )
    if #segList == 0 then 
        -- nothing we can do!
        print("  PushPull failed due to lack of distant segs")
        return false 
    end
    centerSeg = segList[ random( #segList) ]

    print( "Push-Pull for badseg ".. badSeg.." and centerseg "..centerSeg )
    changeSucceeded = PutPushPullBands( badSeg, 5, centerSeg, wantPushNotPull )
    if changeSucceeded then SimpleBandedWiggleWrapper( time, 100.0 ) end
    
    DelBands()
    return changeSucceeded
end

function PickVoidCrusherAction( idx )
    local time = QUICKWIGGLE_ITERS
    DelBands()

    print( "Void crusher for " .. idx )
    changeSucceeded = PutVoidCrusherBands( idx )
    if changeSucceeded then SimpleBandedWiggleWrapper( time, 100.0 ) end

    DelBands()
    return changeSucceeded
end

function PickMultibander( )
    local time = QUICKWIGGLE_ITERS
    changeSucceeded = false

    DelBands()
    local which = random( 3 )
    if which == 1 then
        print("Multiple ZL bands in space")
        changeSucceeded = PutMultipleZLBandsInSpace( 5, 1.0 )
    elseif which == 2 then
        print("Multiple bands between segs")
        changeSucceeded = PutMultipleBandsBetween( 5, 0.5 )
    else -- if which == 3
        print("Multiple bands both types")
        changeSucceeded = PutMultipleBandsBothTypes( 5, 5, 0.5, 1.0 )
    end

    if changeSucceeded then
        -- maybe change to SimpleBandedWiggleWrapper( time, 100.0 ) 
        local which = random( 4 )
        if which == 1 then
            qStabWiggle( time, false )
        elseif which == 2 then
            WiggleWalkNicely( time )
        elseif which == 3 then
            WiggleZlbNicely( time, 1.0 ) -- time, strength
        else -- if which == 4 then
            LocalWiggleNicely( time, 6, 3 ) -- time, chunksize, passcount
        end

        DelBands()

        which = random( 4 )
        if which == 1 then
            qStabWiggle( time, true ) -- big enough change that we *ought* to work harder
        elseif which == 2 then
            WiggleWalkNicely( time )
        elseif which == 3 then
            WiggleZlbNicely( time, 1.0 ) -- time, strength
        else -- if which == 4 then
            LocalWiggleNicely( time, 6, 3 ) -- time, chunksize, passcount
        end
    end

    DelBands()
    return changeSucceeded
end

function PickLSQuaker( idx, wantFreeze )
    if wantFreeze == nil then wantFreeze = true end
    local time = QUICKWIGGLE_ITERS

    if wantFreeze then freeze.UnfreezeAll() end  -- I think we don't want to tangle freezes

    print("PickLSQuaker for ".. idx)
    changeSucceeded =  PutManyBandsToSeg( idx, 0.50, 6, true, true, wantFreeze )
    if changeSucceeded then
        BandedWiggleUntilTargetMet( time, 150.0, 5, CORE_PULLINGCI) -- time, scoredelta, maxtries
        if wantFreeze then freeze.UnfreezeAll() end
        DelBands()
        qStabWiggle( time )
    end
    freeze.UnfreezeAll()  -- shouldn't be necessary, but paranoia is good for the soul
    return changeSucceeded
end

function PickMutatorAlgorithm( idx )
    print( "mutating segment " .. idx )
    r = random()
    if r < PROB_MUTATESIMILAR then
        changeSucceeded = MutateToSimilar(idx)
    else
        changeSucceeded = MutateToAny(idx)
    end
    return changeSucceeded
end

----------------------------------------------------------------------
-- CHOOSERS THAT SORT BY SCALE OF OPERATION
----------------------------------------------------------------------

function PickLocalOperation( idx )
    if idx == nil then idx = randomMovableSeg( ) end
    local time = QUICKWIGGLE_ITERS
    local didIdealize = false
    local glyHingeAllowed = (idx > 1) and (idx < segCt) and (structure.GetAminoAcid( idx ) == 'g')
    local allowMutate = not isNoMutate( idx )

    choiceCt = 3
    if glyHingeAllowed or allowMutate then 
        choiceCt = choiceCt + 1 
    end
 
    local which = random( choiceCt )
    if which == 1 then
        local strength = random( 0.5, 1.5, true )
        print( "BIS for segment " .. idx.." with strength "..strength )
        changeSucceeded = PutBandInSpace( idx, 5.0, strength )  -- 5.0 = maxRho
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, idx, idx, 0.20 ) end
    elseif which == 2 then
        local strength = random( 0.5, 1.5, true )
        print( "Band-Between for segment " .. idx.." with strength "..strength )
        changeSucceeded = PutBandToRandomSeg( idx, strength )
        if changeSucceeded then SimpleWiggleBeforeAndAfter( time, idx, idx, 0.20 ) end
    elseif which == 3 then
        if USE_CUTS then
            print( "IdealizeInSmallRegion with cuts for "..idx )
        else
            print( "IdealizeInSmallRegion withOUT cuts for "..idx )
        end
        changeSucceeded = IdealizeInSmallRegion( idx, time, USE_CUTS )
        if changeSucceeded then qStabWiggle( time, random() < 0.50 ) end
    else -- if which == 4 then
        if glyHingeAllowed and allowMutate then
             which = random(2)
             if which == 1 then glyHingeAllowed = false end
        end

        if glyHingeAllowed then
            print( "Glycine-hinge for segment " .. idx )
            changeSucceeded = PutGlycineHingeBands( idx, random( 0.25, 0.75 ) ) -- bandstrength
            if changeSucceeded then SimpleWiggleBeforeAndAfter( time, idx, idx, 0.0 ) end
        else 
            print( "Spot-mutation for segment " .. idx )
            changeSucceeded = PickMutatorAlgorithm( idx ) 
            if changeSucceeded then SimpleWiggleBeforeAndAfter( time, idx, idx, 0.90 ) end
        end
    end

    freeze.UnfreezeAll()  -- if we froze anything, we should clean that up now
    DelBands()  -- if we put in bands, they should go now
    if not changeSucceeded then print( "CHANGE FAILED" ) end
    return changeSucceeded
end

function PickMidsizeOperation( )
    changeSucceeded = false

    local idx, startIdx, endIdx 

    local which = random( 6 )
    if which == 1 then
        changeSucceeded = PickMultibander( )
    elseif which == 2 then
        idx = randomMovableSeg()
        local startIdx, endIdx = getRegionForOperation( idx )
        if endIdx - startIdx > 1 then
            changeSucceeded = PickRegionBander( startIdx, endIdx )
        end
    elseif which == 3 then
        idx = randomMovableSeg()
        changeSucceeded = PickPushPullAction( idx, random() < 0.50 )
     elseif which == 4 then
        idx = randomMovableSeg()
        changeSucceeded = PickVoidCrusherAction( idx )
    elseif which == 5 then
        idx = randomMovableSeg()
        changeSucceeded = PickLSQuaker( idx, true ) -- or randomly vary true/false
    else -- if which == 6 then
       idx = randomMovableSeg()
       changeSucceeded = PickPointRebuildOperation( idx, 5 )   -- these are pricy, limit to 5 attempts
    end

    freeze.UnfreezeAll()  -- if we froze anything, we should clean that up now
    DelBands()  -- if we put in bands, they should go now
    if not changeSucceeded then print( " CHANGE FAILED" ) end
    return changeSucceeded
end

function PickLargeScaleOperation( )
    local time = QUICKWIGGLE_ITERS
    local maxStren = 0.75

    local which = random(2)
    if which == 1 then
        changeSucceeded = PickCentralBander( random(5), maxStren ) -- maxCores
    else -- if which == 2 then
        changeSucceeded = PickEveryXToEveryYBanding( maxStren ) 
    end

    freeze.UnfreezeAll()  -- if we froze anything, we should clean that up now
    DelBands()  -- if we put in bands, they should go now
    if not changeSucceeded then print( " CHANGE FAILED" ) end
    return changeSucceeded
end

function PickTargetedOperation( idx, probLocal )
    local time = QUICKWIGGLE_ITERS
    local didIdealize = false

    if idx == 0 then
        idx = randomLowScoringSeg( random(10) )
    end
    if probLocal == nil then 
        probLocal = PROB_LOCAL / ( PROB_MIDSIZE + PROB_LOCAL )
    end

    local glyHingeAllowed = (idx > 1) and (idx < segCt) and (structure.GetAminoAcid( idx ) == 'g')
    local allowMutate = not isNoMutate( idx )

    local r = random()
    if r < probLocal then
        changeSucceeded = PickLocalOperation( idx )
     else
        -- do a targetable midsize action (ie leave out multibanders)
        local startIdx, endIdx
        local which = random( 5 )

        if which == 1 then
            startIdx, endIdx = getRegionForOperation( idx )
            if endIdx - startIdx > 1 then  -- don't really expect to ever fail, but...
                changeSucceeded = PickRegionBander( startIdx, endIdx )
            end
        elseif which == 2 then
            changeSucceeded = PickPushPullAction( idx, random() < 0.50 )
        elseif which == 3 then
            changeSucceeded = PickVoidCrusherAction( idx )
        elseif which == 4 then
            changeSucceeded = PickLSQuaker( idx, true ) -- or randomly vary true/false
        else -- if which == 5 then
            changeSucceeded = PickPointRebuildOperation( idx, 5 )   -- these are pricy, limit to 5 attempts
        end
    end
    return changeSucceeded
end

--------------------------------------------------------------------------------------------------------
--
---- THE TOP-LAYER PROTEIN MODIFIERS: 
--
--------------------------------------------------------------------------------------------------------

function DoCleanupWiggle( time, minGainForIteration, rbscore, qsHolder )
    -- look for any transient high score, grab it
    local rbnewscore = getScore( true )
    sawTransient = false
    if rbscore ~= rbnewscore and rbnewscore > getScore() then
        print( "    transient score: "..rbnewscore)
        sawTransient = true
        save.Quicksave( qsHolder )   -- this is the "produced neighbor"
        recentbest.Restore()              -- "current" has the transient

        -- transients tend to need cleaning, but delay removing cuts till *after* re-wiggle
        DelBands()
        selection.DeselectAll()
        freeze.UnfreezeAll()
    end

     -- Use "slow" wiggle functions
    if sawTransient then
        -- try very hard to get a good score
        FuseWiggleNicely( time )   -- expensive time-wise; look for better choices
        for i=1, segCt do
            RemoveCutAtIndex( i ) -- we want to be rid of cuts, even at a point loss
        end

        if getScore() < rbnewscore then
            -- TODO consider trying to track scores more closely to see if cleaning a cut is why 
            -- we lost points. If that is the case, then try to fight to regain them (rebuilds, idealize)
            save.Quickload( qsHolder )
            sawTransient = false   -- give ourselves a second chance!
        end
    end

    if not sawTransient then
        ShakeOrMutate( time/2 )

        local r = random( ) 
        if r < PROB_SLOWWIGGLE_ALL then
            print( "IterativeWiggleAll")
            IterativeWiggleAll( time, minGainForIteration, true, true) -- a global wiggler
        elseif r < PROB_SLOWWIGGLE_ALL + PROB_SLOWWIGGLE_LOCAL then
            print( "IterativeWiggleLocalByChunk" )  -- can be very slow if there's a big score change to fix!
            IterativeWiggleLocalByChunk( time, 1, segCt, minGainForIteration, random(3,6), random()<0.50 )
        else
            print( "IterativeCutAndWiggle")   -- slow enough that we want to do this more rarely
            IterativeCutAndWiggle( time, minGainForIteration, random(3,13) )
        end

        for i=1, segCt do              -- not *supposed* to be necessary here, but
            RemoveCutAtIndex( i ) -- we want to be sure to be rid of any stray cuts.
        end
    end
end

function IsNeighborDifferentFromOld( oldScoreMat, threshold )
    local newScoreMat = {}
    local scoreVector = {}
    ConstructScoreMatrix( newScoreMat )
    GetEuclideanDistanceVector( oldScoreMat, newScoreMat, scoreVector )
    for i=1, #scoreVector do
        if scoreVector[i] > threshold then return true end
    end
    return false
end

function IsProteinChanged( oldScore, oldScoreMat, threshold )
    if not current.AreConditionsMet() then return false end
    local currScore = getScore( )
    if currScore > oldScore then return true end
    return IsNeighborDifferentFromOld( oldScoreMat, threshold )
end

function GetTargetedNeighborState(idx, threshold, wantCleanupWiggle)
    local oldScoreMat = {}
    local oldScore = getScore( )
    ConstructScoreMatrix( oldScoreMat )

    if wantCleanupWiggle  == nil then wantCleanupWiggle  = true end
    DelBands()
    selection.DeselectAll()
    freeze.UnfreezeAll()

    recentbest.Save()
    local rbscore = getScore( true )

    print( "Targeting "..idx)
    save.Quicksave( QS_NeighborTemp )
    local goodNeighbor = false
    while not goodNeighbor do
        local changeSucceeded = false
        goodNeighbor = false

        save.Quickload( QS_NeighborTemp )

        changeSucceeded = PickTargetedOperation( idx )
        if changeSucceeded then
            goodNeighbor = IsProteinChanged( oldScore, oldScoreMat, threshold )
            if not goodNeighbor then
                print( "TARGNEIGHBOR REJECTED - redoing" )
            end
        end
    end   -- END while not goodNeighbor

    if wantCleanupWiggle then
        local time = SLOWWIGGLE_ITERS
        DoCleanupWiggle( time, SLOW_IT_WIGGLE_GAIN, rbscore, QS_NeighborTemp )
    end

    print( "TargetedNeighbor Score: " .. getScore())
end

function GetNeighborState( threshold, wantCleanupWiggle )
    local oldScoreMat = {}
    local oldScore = getScore( )
    ConstructScoreMatrix( oldScoreMat )

    if wantCleanupWiggle  == nil then wantCleanupWiggle  = true end
    DelBands()
    selection.DeselectAll()
    freeze.UnfreezeAll()

    recentbest.Save()
    local rbscore = getScore( true )

    save.Quicksave( QS_NeighborTemp )
    local goodNeighbor = false
    while not goodNeighbor do
        local changeSucceeded = false
        goodNeighbor = false

        save.Quickload( QS_NeighborTemp )
        local r = random( )   -- 0.25 appears too high. try #2 = 0.10. If still too high, then drop to 0.05. If still, then abandon targeting
        if r < PROB_TARGETED then
            changeSucceeded = PickTargetedOperation( 0 )
        elseif r < PROB_TARGETED + PROB_LOCAL then
            changeSucceeded = PickLocalOperation()
        elseif r < PROB_TARGETED + PROB_LOCAL + PROB_MIDSIZE then
            changeSucceeded = PickMidsizeOperation()
        else 
            changeSucceeded = PickLargeScaleOperation()
        end

        if changeSucceeded then
            goodNeighbor = IsProteinChanged( oldScore, oldScoreMat, threshold )
            if not goodNeighbor then
                print( "NEIGHBOR REJECTED (score:"..getScore()..")")
            end
        end
    end   -- END while not goodNeighbor

    if wantCleanupWiggle then
        local time = SLOWWIGGLE_ITERS
        DoCleanupWiggle( time, SLOW_IT_WIGGLE_GAIN, rbscore, QS_NeighborTemp )
    end

    DebugPrint( "Neighbor Score: " .. getScore())
end

----------------------------------------------------------------------
--
---- SIMULATED ANNEALING ALGORITHM FUNCTIONS
--
----------------------------------------------------------------------
function temperature(num) 
    --      num goes 0 -> 1, so this goes +INF -> 0 as 1/x
    -- Real research could be done (rosetta folks have done) on what a "good"
    -- temperature function is. 1/x seems to do well enough for me.
    return ( 1.0 / (num + 0.000001) ) - 1.0
end

function P(ecurr, enew, T)  -- this is standard in SA work
    if enew < ecurr then return 1.0 end
    return math.exp( (ecurr - enew)/T )
end

function E()  -- must negate score to function in SA. 8000 is irrelevant, just for scaling
    return 8000.0 - getScore()
end

function SA( kmax, emax, failsAllowed, qs_slot )  -- emax is expected to "never happen"
    local ecurr = E()
    save.Quicksave( QS_SATemp )
    local ebest = ecurr 

    local k = 0.0 
    local failcount = 0
    while k < kmax and ecurr > emax do
        save.Quickload( QS_SATemp )

        local T = temperature(k/kmax)
        print("##"..k..": temp="..TrimNum(T).." fails="..failcount.." score="..TrimNum( getScore()) )
        GetNeighborState( 0.1, true )  -- TODO: fool with threshold. (so far, seems good)
        local enew = E()
        local p = P(ecurr, enew, T)

        local r = random()
        if p > r then
            print( ">>>> accepting state with score " .. TrimNum( getScore()) )
            save.Quicksave( QS_SATemp )
            ecurr = enew 
            failcount = 0
        else
            failcount = failcount + 1
            if failcount > failsAllowed then   
                return  -- give up on this run
            end
        end
        if enew < ebest then 
            save.Quicksave( qs_slot )
            ebest = enew 
        end
        k = k + 1.0 
    end
end

----------------------------------------------------------------------
-- ALGORITHM OPTIONS AVAILABLE FROM MAIN
----------------------------------------------------------------------

function PerformSA( qs_start, qs_run, qs_best, preferReheat, failsAllowed, ctRestarts, stepsPerRun, minGainPerRun )
    if minGainPerRun == nil then minGainPerRun = 0 end
    local startScore = getScore()
    local bestScore = startScore
    local startTime = os.time()
    local newscore = 0.0

    for i=1, ctRestarts do
        if preferReheat then
            save.Quickload( qs_run )
        else
            save.Quickload( qs_start )
        end
        local currscore = getScore()
        local currtime = os.time()
        print( "START: SA#"..i.." score: "..TrimNum(currscore) )

        SA( stepsPerRun, -999999.0, failsAllowed, qs_run )

        save.Quickload( qs_run )
        newscore = getScore()
        gain = newscore - bestScore
        if( newscore > bestScore ) then
            save.Quicksave( qs_best )
            bestScore = newscore
        end
        print( "FINISH: SA#"..i.." score: "..TrimNum(newscore)..
               " gain: "..newscore-currscore.." time="..os.time() - currtime )
        if gain < minGainPerRun then
            print("  Insufficient Gain "..gain.." - halting SA")
            break
        end
    end
    print( "SA: steps: ".. ctRestarts*stepsPerRun ..
             " gain: "..newscore-startScore.."  time: "..os.time()-startTime)
end

function MarchForward( qs_best, maxRunsNoGain, stepsWanted, neighborThreshold )
    if neighborThreshold == nil then neighborThreshold = 0.01 end

    local ctRunsWithoutGain = 0    
    local startScore = getScore()
    local bestScore = getScore()
 
    recentbest.Save()
    for i=1, stepsWanted do
         if ctRunsWithoutGain >= maxRunsNoGain then return end

        save.Quickload(qs_best) 
        local prevscore = getScore()
        print( "MF STEP ".. i.. " current score: ".. prevscore)

        GetNeighborState( neighborThreshold, true )
        local currscore = getScore()
        if (bestScore < currscore) then
            ctRunsWithoutGain = 0
            gain = currscore - bestScore
            bestScore = currscore
            save.Quicksave(qs_best)
            print( ">>>> Gain of ".. gain .. " points" )
        else
            gain = 0
            ctRunsWithoutGain=ctRunsWithoutGain+1
            recentbest.Restore()
        end
    end
    print( "MF FINISH: total gain ="..bestScore-startScore )
end

-- consider supporting 3 modes: random idx, count from 1 to segCt, and "by worst"
-- ("by worst" is currently the only supported mode) 
-- AS IT STANDS, THIS IS WAY TOO MUCH TIME FOR WAY TOO LITTLE GAIN. It might 
-- benefit by attacking neighborhoods of low-scoring areas or by adding targeted wiggles or by bouncing
-- among various low-scoring choices...
function TargetedMarchForward( qs_best, maxRunsNoGain, maxStepsPerAttempt, neighborThreshold )
    if neighborThreshold == nil then neighborThreshold = 0.01 end

    probLocal = PROB_LOCAL / (PROB_LOCAL + PROB_MIDSIZE)
    local ctRunsWithoutGain = 0    
    local startScore = getScore()
    local bestScore = getScore()
    local notDoneYet = true

    local allSegs = {}
    local segsVisited = {}
    for i=1, segCt do
        allSegs[ i ] = i
        segsVisited[ i ] = not isMovableSeg( i )  -- no point visiting non-movable segments!
    end

    while notDoneYet do
        -- pick an index:
        getWorstScoringAAs( allSegs, segCt )
        k = 1
        while segsVisited[ allSegs[ k ] ] == true and k <= segCt do k = k + 1 end
        if k == segCt + 1 then break end  -- nothing left!
        segsVisited[ allSegs[ k ] ] = true
        idx = allSegs[ k ]

        -- begin the march!
        print("  Targeting "..idx.." for "..maxStepsPerAttempt.." tries")
        recentbest.Save()
        ctRunsWithoutGain = 0
        for i=1, maxStepsPerAttempt do
             if ctRunsWithoutGain >= maxRunsNoGain then break end

            save.Quickload(qs_best) 
            local prevscore = getScore()
            print( "MF-T STEP ".. i.. " current score: ".. prevscore)

            GetTargetedNeighborState( idx, probLocal, neighborThreshold, true )

            local currscore = getScore()
            if (bestScore < currscore) then
                ctRunsWithoutGain = 0
                gain = currscore - bestScore
                bestScore = currscore
                save.Quicksave(qs_best)
                print( ">>>> Gain of ".. gain .. " points" )
            else
                gain = 0
                ctRunsWithoutGain=ctRunsWithoutGain+1
                recentbest.Restore()
            end
        end
    end
end

function PerformMarchWiggleAction( qs_start, qs_run, qs_best, lastScore, marchThreshold, whichWiggle )
    save.Quickload( qs_best )
    newscore = getScore()
    print( "MarchWiggle: lastscore="..lastScore..", currscore="..newscore )

    -- see if we should attempt some MarchForward
    if newscore - lastScore > marchThreshold then
        PerformMajorOperation( qs_start, qs_run, qs_best, "MarchForward" )
        lastScore = getScore()
    end

    print("lastScore = "..lastScore)

    -- perform a default version of whatever wiggle the user asked for
    local time = FINALWIGGLE_ITERS
    local gain = FINALRUN_WIGGLE_GAIN
    print( "   "..whichWiggle )
    if whichWiggle == "IterativeWiggleLocalByChunk-freeze" then
        IterativeWiggleLocalByChunk( time, 1, segCt, gain, random(3,5), true )
    elseif whichWiggle == "IterativeWiggleLocalByChunk" then
        IterativeWiggleLocalByChunk( time, 1, segCt, gain, random(3,5), false )
    elseif whichWiggle == "FreezeBurst" then
        PerformFreezeBurst( time )
    elseif whichWiggle == "IterativeWiggleAll" then
        IterativeWiggleAll( time, gain, true, true)
    elseif whichWiggle == "MonsterFuse" then
        FuseWiggleNicely( time, true )
    elseif whichWiggle == "FuseBurstIterateCi" then
        FuseBurstIterateCi ( time )
    elseif whichWiggle == "FullIterativeCutAndWiggle" then
        FullIterativeCutAndWiggle( time, gain, 17, true)
    else
        print("!! no such wiggle")
    end

    if getScore() > newscore then 
        DelBands( )  -- shouldn't be any, but just to be safe...
        freeze.UnfreezeAll( )
        selection.DeselectAll( )
        save.Quicksave( qs_best ) 
    else
         save.Quickload( qs_best )
    end

    return lastScore
end

function PerformMarchWiggleSequence( qs_start, qs_run, qs_best, algSequence )
    print( " Begin MarchWiggle sequence" )
    if algSequence == nil or #algSequence == 0 then 
        algSequence = {"IterativeWiggleLocalByChunk-freeze", "FreezeBurst", 
                                "IterativeWiggleAll", "MonsterFuse", 
                                "IterativeWiggleLocalByChunk", 
--  FuseBurstIterateCi disabled by LociOiling 2014/04/27
--                                "FuseBurstIterateCi", 
                                "FullIterativeCutAndWiggle"}
    end

    local lastScore = getScore()   -- first call won't do march forward
    local minGain = MARCHWIGGLE_MINGAIN
    for i=1, #algSequence do
        lastScore = PerformMarchWiggleAction( qs_start, qs_run, qs_best, lastScore, minGain, algSequence[ i ] )
    end
end

function PerformDefaultSequence( qs_start, qs_run, qs_best )
    save.Quicksave( qs_best )
    save.Quicksave( qs_start )
    save.Quicksave( qs_run )

    -- start with a classic Simulated Annealing sequence
    PerformMajorOperation( qs_start, qs_run, qs_best, "SA" )

    PerformMarchWiggleSequence( qs_start, qs_run, qs_best )
end

function PerformMajorOperation( qs_start, qs_run, qs_best, whichOperation )
    print( "Operation Start: Score="..getScore() )

    save.Quicksave( qs_start )
    save.Quicksave( qs_run)
    save.Quicksave( qs_best )

    if whichOperation == "SA" then
        print( "SA: Using " .. SA_STEPSPERRUN .. " steps for " .. SA_NUMSTARTS .. " runs" )
        PerformSA( qs_start, qs_run, qs_best, 
                        REHEAT_NOT_RESTART, SA_FAILSBEFORERESTART, 
                        SA_NUMSTARTS, SA_STEPSPERRUN, SA_MINGAIN)
    elseif whichOperation == "DefaultSequence" then
        print( "DefaultSeqence: Score="..getScore() )
        PerformDefaultSequence( qs_start, qs_run, qs_best )
    elseif whichOperation == "MarchForward" then
        print( "MF: Using " .. MF_STEPSPERRUN .. " steps. Score="..getScore() )
        MarchForward( qs_best, MF_FAILSBEFORERESTART, MF_STEPSPERRUN ) 
    elseif whichOperation == "TargetedMarchForward" then
        print( "MF-T: Using " .. TARG_STEPSPERRUN .. " steps. Score="..getScore() )
        TargetedMarchForward( qs_best, TARG_FAILSBEFORERESTART, TARG_STEPSPERRUN )
    elseif whichOperation == "MarchWiggle" then
        print( "MarchWiggle" )
        PerformMarchWiggleSequence( qs_start, qs_run, qs_best )  -- support an algSequence here!
    else
        print( "!! Unknown operation" )
    end

    save.Quickload( qs_best )
    print( "Operation Finished: Score="..getScore() )
end

--------------------------------------------------------------------------
-- SETUP AND CLEANING
-- TODO: Improve dialog to allow user to tweak various times, strengths,
-- algorithm choices (and to set "don't mutate idx" choices), and more
--------------------------------------------------------------------------
function cleanup()
    --
    -- restore original protein state
    --
    behavior.SetClashImportance( InitialClashImportance )
    selection.DeselectAll()
    freeze.UnfreezeAll()
    for i=1,segCt do
        RemoveCutAtIndex( i ) -- we can try, at least!
    end
    for i=1,segCt do
        structure.SetSecondaryStructure( i, InitialSecondaryIdxes[i] )
    end
    for i=1,segCt do
        if InitialSelectedIdxes[i] then
            selection.Select( i )
        end
    end
    for i=1,segCt do
        bFroz, sFroz = InitialFrozenTable[i]
        if bFroz or sFroz then
            freeze.Freeze( i, bFroz, sFroz)
        end
    end

    -- do what I can with bands. Cannot simply re-enable original bands
    -- because can't have them sitting around disabled while working.
    -- Someday try disabling all bis bands (then move protein?) and recording
    -- all between-seg bands for later restoration. (if goal length != 3.50, then
    -- set goal length)
    DelBands()
    --InitialBandCount = band.GetCount()  -- restore bands
    --for i=1, InitialBandCount do
    --    if InitialEnabledBandTable[i] then band.Enable( i )
    --end
end

function PrintState( )
    local gain = getScore() - InitialScore
    if gain < 0.001 then 
        print( " No change" )
    else
        print( " Startscore: "..TrimNum( InitialScore ))
        print( " Score: "..TrimNum( getScore() ) )
        print( " Total gain: "..TrimNum( gain ))
    end
    print("")
    print( " Run time "..os.time() - StartTime)
end

function End(errstr)
    if EndCalled then return end    -- no infinite recursion please!!!
    EndCalled = true

    print("")
    if errstr ~= nil then
        if string.find(errstr, "Cancelled") then
            print( "User cancel" )
        else
            print( errstr )
        end
    end

    save.Quickload( QS_Best )
    PrintState( )

    cleanup()
end

function InitializePuzzleState()
    seedRandom()

    -- learn a few things
    FindAllMovableSegs()
    getHelices( )
    getSheets( )
    getLoops()
    PuzzleAllowsMutate = IsPuzzleMutable( )
    if PuzzleAllowsMutate then IsDesignerPuzzle = checkAla() end
    InitialScore = getScore( )

    --
    -- store original protein state and set to "clean" state
    --
    InitialClashImportance = behavior.GetClashImportance()
    behavior.SetClashImportance(1.0)

    for i=1, segCt do     -- save "secondary structures", selectedness, frozennesses
        InitialSecondaryIdxes[i] = structure.GetSecondaryStructure( i )
        InitialSelectedIdxes[i] = selection.IsSelected(i)
        local bFroz, sFroz = freeze.IsFrozen(i)
        InitialFrozenTable[i] = bFroz, sFroz
    end
    selection.DeselectAll()
    freeze.UnfreezeAll()

    --InitialBandCount = band.GetCount()   --save bands
    --for i=1, InitialBandCount do
    --    InitialEnabledBandTable[i] = band.IsEnabled( i )
    --    band.Disable( i )
    --end
    DelBands()
    recentbest.Save()
end

-- A DIALOG BOX SEQUENCE FROM THE NETHERWORLD
function LetUserChooseParams( )
    dlg1 = dialog.CreateDialog("Simulated Annealing And More")
    dlg1.cmt1 = dialog.AddLabel("Pick one of following (topmost checked will be run)")
    dlg1.wantDef = dialog.AddCheckbox("Perform Default Sequence: (SA, then March-Wiggles)", true)
    dlg1.wantJustSA = dialog.AddCheckbox("Perform only SA", false)
    dlg1.wantJustMW = dialog.AddCheckbox("Perform only March-Wiggles", false)
--    dlg1.wantJustWig = dialog.AddCheckbox("Perform only Wiggles", false)
--    dlg1.wantJustBander = dialog.AddCheckbox("Perform simple actions", false)
    dlg1.ok = dialog.AddButton("OK", 1)
    dlg1.cancel = dialog.AddButton("Cancel", 0)
    if dialog.Show(dlg1) == 0 then return false end

    if dlg1.wantDef.value then whichAlg = "DefaultSequence"
    elseif dlg1.wantJustSA.value then whichAlg = "SA"
    elseif dlg1.wantJustMW.value then whichAlg = "MarchWiggle"
    elseif dlg1.wantJustWig.value then whichAlg = "Wiggle"
    elseif dlg1.wantJustBander.value then whichAlg = "Changer"
    else return false -- something went wrong...
    end

    if whichAlg == "DefaultSequence" or whichAlg == "SA" then
        dlg1 = dialog.CreateDialog("Simulated Annealing Parameters")
        dlg1.stepsper = dialog.AddSlider( "StepsPerRun", SA_STEPSPERRUN, 25, 200, 0 )
        dlg1.starts = dialog.AddSlider( "RunStarts", SA_NUMSTARTS, 1, 10, 0 ) 
        dlg1.samingain = dialog.AddSlider( "MinGainForRun", SA_MINGAIN, 0, 5, 1 )
        dlg1.advancedpars = dialog.AddCheckbox("Advanced parameters", false)
        dlg1.ok = dialog.AddButton("OK", 1)
        dlg1.cancel = dialog.AddButton("Cancel", 0)
        if dialog.Show(dlg1) == 0 then return false end

        SA_STEPSPERRUN = dlg1.stepsper.value
        SA_NUMSTARTS = dlg1.starts.value
        SA_MINGAIN = dlg1.samingain.value
        if dlg1.advancedpars.value then 
            dlg1 = dialog.CreateDialog("Advanced Simulated Annealing Parameters")
            dlg1.sep2 = dialog.AddLabel( "----------------Scales of Action of Steps---------------------" )
            dlg1.comment1 = dialog.AddLabel( "Scale actions are 'probabilities' - scaled to sum to 1" )
            dlg1.scaleLocal = dialog.AddSlider( "S-Local", PROB_LOCAL, 0.01, 1, 2 )
            dlg1.scaleMid = dialog.AddSlider( "S-Midsize", PROB_MIDSIZE, 0.01, 1, 2 )
            dlg1.scaleTarg = dialog.AddSlider( "S-Targeted", PROB_TARGETED, 0, 1, 2 )
            dlg1.scaleLarge= dialog.AddSlider( "S-Largescale", 
                            1.0-(PROB_LOCAL+PROB_MIDSIZE+PROB_TARGETED), 0, 1, 2 )
            dlg1.sep3 = dialog.AddLabel( "----------------Quick Wiggle During Step---------------------" )
            dlg1.qwtime = dialog.AddSlider( "QuickWiggleIters", QUICKWIGGLE_ITERS, 2, 10, 0)
            dlg1.qwgain = dialog.AddSlider( "QuickWiggleGain", QUICK_IT_WIGGLE_GAIN, 0, 1, 2)
            dlg1.sep4 = dialog.AddLabel( "----------------Slow Wiggle After Step------------------------" )
            dlg1.swtime = dialog.AddSlider( "SlowWiggleIters", SLOWWIGGLE_ITERS, 2, 14, 0 )
            dlg1.swgain = dialog.AddSlider( "SlowWiggleGain", SLOW_IT_WIGGLE_GAIN, 0, 1, 2)
            dlg1.comment3 = dialog.AddLabel( "Wiggle choices are 'probabilities' - scaled to sum to 1" )
            dlg1.wall = dialog.AddSlider( "W-Choice All", PROB_SLOWWIGGLE_ALL, 0.01, 1, 2 )
            dlg1.wloc = dialog.AddSlider( "W-Choice Local", PROB_SLOWWIGGLE_LOCAL, 0.01, 1, 2 )
            dlg1.wfuse= dialog.AddSlider( "W-Choice Fuse", 
                               1.0-(PROB_SLOWWIGGLE_ALL+PROB_SLOWWIGGLE_LOCAL), 0, 1, 2 )
            dlg1.usecuts = dialog.AddCheckbox("Allow cuts during idealize", USE_CUTS)
            dlg1.musthealcuts = dialog.AddCheckbox("Forbid unhealed cuts after idealize", FORBID_UNHEALED_CUTS)

            if whichAlg == "DefaultSequence" then
                dlg1.sep5 = dialog.AddLabel( "-----------------After-Run March-Wiggles-----------------------" )
                dlg1.mfstepsper = dialog.AddSlider( "StepsPerRun", MF_STEPSPERRUN, 5, 100, 0 )
                dlg1.mffails = dialog.AddSlider( "FailsBeforeQuit", MF_FAILSBEFORERESTART, 5, 100, 0 ) 
                dlg1.fwtime = dialog.AddSlider( "Wiggle Iters", FINALWIGGLE_ITERS, 2, 20, 0 )
                dlg1.comment2 = dialog.AddLabel( "WiggleGain is decimal places: 1, 0.1, .., 0.0001")
                dlg1.fwgain = dialog.AddSlider( "Wiggle Gain", 3, 0, 4, 0 ) -- there is no good constant to use here!
                dlg1.mwmingain = dialog.AddSlider( "Remarch Gain", MARCHWIGGLE_MINGAIN, 0, 5, 1 )
            end
            dlg1.ok = dialog.AddButton("OK", 1)
            dlg1.cancel = dialog.AddButton("Cancel", 0)
 
            if dialog.Show(dlg1) == 0 then return false end

            QUICKWIGGLE_ITERS = dlg1.qwtime.value
            SLOWWIGGLE_ITERS = dlg1.swtime.value
            QUICK_IT_WIGGLE_GAIN = dlg1.qwgain.value
            SLOW_IT_WIGGLE_GAIN = dlg1.swgain.value
            local t1 = dlg1.scaleLocal.value
            local t2 = dlg1.scaleMid.value
            local t3 = dlg1.scaleTarg.value
            local t4 = dlg1.scaleLarge.value
            PROB_LOCAL = t1 / (t1 + t2 + t3 + t4)
            PROB_MIDSIZE = t2 / (t1 + t2 + t3 + t4)
            PROB_TARGETED = t3 / (t1 + t2 + t3 + t4)
            t1 = dlg1.wall.value
            t2 = dlg1.wloc.value
            t3 = dlg1.wfuse.value
            PROB_SLOWWIGGLE_ALL = t1 / (t1 + t2 + t3)
            PROB_SLOWWIGGLE_LOCAL = t2 / (t1 + t2 + t3)
            USE_CUTS = dlg1.usecuts.value
            FORBID_UNHEALED_CUTS = dlg1.musthealcuts.value
            DebugPrint("QuickIters="..QUICKWIGGLE_ITERS)
            DebugPrint("SlowIters="..SLOWWIGGLE_ITERS )
            DebugPrint("QuickGain="..QUICK_IT_WIGGLE_GAIN)
            DebugPrint("SlowGain="..SLOW_IT_WIGGLE_GAIN)
            DebugPrint("Prob_local="..PROB_LOCAL)
            DebugPrint("Prob_midsize="..PROB_MIDSIZE)
            DebugPrint("Prob_targetd="..PROB_TARGETED)
            DebugPrint("Prob_slowwiggle-all="..PROB_SLOWWIGGLE_ALL)
            DebugPrint("Prob_slowwiggle-local="..PROB_SLOWWIGGLE_LOCAL)
            if USE_CUTS then DebugPrint(" cuts are allowed") else DebugPrint(" cuts are forbidden") end
            if FORBID_UNHEALED_CUTS then DebugPrint(" cuts must get healed") else DebugPrint(" cuts can stay unhealed") end
            if whichAlg == "DefaultSequence" then
                MF_STEPSPERRUN = dlg1.mfstepsper.value
                MF_FAILSBEFORERESTART = dlg1.mffails.value
                FINALWIGGLE_ITERS = dlg1.fwtime.value
                MARCHWIGGLE_MINGAIN = dlg1.mwmingain.value
                if dlg1.fwgain.value == 0 then 
                    FINALRUN_WIGGLE_GAIN = 1 
                elseif dlg1.fwgain.value == 1 then
                    FINALRUN_WIGGLE_GAIN= 0.1
                elseif dlg1.fwgain.value == 2 then
                    FINALRUN_WIGGLE_GAIN = 0.01
                elseif dlg1.fwgain.value == 3 then
                    FINALRUN_WIGGLE_GAIN = 0.001
                elseif dlg1.fwgain.value == 4 then
                    FINALRUN_WIGGLE_GAIN = 0.0001
                end
                DebugPrint("MarchForward steps="..MF_STEPSPERRUN)
                DebugPrint("MarchForward failsallowed="..MF_FAILSBEFORERESTART)
                DebugPrint("FinalWiggle iters="..FINALWIGGLE_ITERS)
                DebugPrint("MarchWiggle mingain="..MARCHWIGGLE_MINGAIN)
                DebugPrint("FinalWiggle mingain="..FINALRUN_WIGGLE_GAIN)
            end
        end -- IF "advance parameters for SA/MW chosen"
    elseif whichAlg == "MarchWiggle" then
        dlg1 = dialog.CreateDialog("March-Wiggle Parameters")
        dlg1.sep5 = dialog.AddLabel( "-----------------------March-Wiggles--------------------------" )
        dlg1.mfstepsper = dialog.AddSlider( "StepsPerRun", MF_STEPSPERRUN, 5, 100, 0 )
        dlg1.mffails = dialog.AddSlider( "FailsBeforeQuit", MF_FAILSBEFORERESTART, 5, 100, 0 ) 
        dlg1.fwtime = dialog.AddSlider( "Wiggle Iters", FINALWIGGLE_ITERS, 2, 20, 0 )
        dlg1.comment2 = dialog.AddLabel( "Wiggle Gain is decimal places: 1, 0.1, .., 0.0001")
        dlg1.fwgain = dialog.AddSlider( "Wiggle Gain", 3, 0, 4, 0 ) -- there is no good constant to use here!
        dlg1.mwmingain = dialog.AddSlider( "Remarch Gain", MARCHWIGGLE_MINGAIN, 0, 5, 1 )
        dlg1.ok = dialog.AddButton("OK", 1)
        dlg1.cancel = dialog.AddButton("Cancel", 0)
        if dialog.Show(dlg1) == 0 then return false end

        MF_STEPSPERRUN = dlg1.mfstepsper.value
        MF_FAILSBEFORERESTART = dlg1.mffails.value
        FINALWIGGLE_ITERS = dlg1.fwtime.value
        MARCHWIGGLE_MINGAIN = dlg1.mwmingain.value
        if dlg1.fwgain.value == 0 then 
            FINALRUN_WIGGLE_GAIN = 1 
        elseif dlg1.fwgain.value == 1 then
            FINALRUN_WIGGLE_GAIN= 0.1
        elseif dlg1.fwgain.value == 2 then
            FINALRUN_WIGGLE_GAIN = 0.01
        elseif dlg1.fwgain.value == 3 then
            FINALRUN_WIGGLE_GAIN = 0.001
        elseif dlg1.fwgain.value == 4 then
            FINALRUN_WIGGLE_GAIN = 0.0001
        end
        DebugPrint("MarchForward steps="..MF_STEPSPERRUN)
        DebugPrint("MarchForward failsallowed="..MF_FAILSBEFORERESTART)
        DebugPrint("FinalWiggle iters="..FINALWIGGLE_ITERS)
        DebugPrint("MarchWiggle mingain="..MARCHWIGGLE_MINGAIN)
        DebugPrint("FinalWiggle mingain="..FINALRUN_WIGGLE_GAIN)
--    elseif whichAlg == "Wiggle" then
--    elseif whichAlg == "Changer" then
    end
    return true
end

function isMajorOperation( whichAlg )
    if whichAlg == "DefaultSequence" then return true end
    if whichAlg =="SA" then return true end 
    if whichAlg == "MarchForward" then return true end
    if whichAlg == "TargetedMarchForward" then return true end
    if whichAlg == "MarchWiggle" then return true end
    return false
end

function main()
    -- some initialization
    InitializePuzzleState()
    didSomething = false

    ------------------------DEBUG STUFF
    -- if any actions could stand some debugging, put them here
    -- then multi-line comment from below to end of function
    ------------------------

    while LetUserChooseParams( ) do
         if isMajorOperation( whichAlg) then
             PerformMajorOperation( QS_Start, QS_Run, QS_Best, whichAlg )
             didSomething = true
--        elseif whichAlg == "Modifier" then
--        elseif whichAlg == "WiggleOnly" then
        end
        save.Quickload( QS_Best )
        PrintState( )
    end

    print( " Done!" )
    if didSomething then End( ) end
end

xpcall( main, End )

Comments

[LociOiling]

This version removes the attempt to run the "FuseBurstIterateCi" function, which was missing in v1.0, resulting in an eventual crash.

KarenCH has plans for a new version of Simulated Annealing, stay tuned….