Recipe: Contact Mapper v2.1

created by KarenCH

Profile

• Name: Contact Mapper v2.1
• ID: 49712
• Shared with: Public
• Parent: None
• Children:
  • Contact Mapper v2.1 no un freeze
  • Contact Mapper v2.1 keep freeze&bands
• Created on: July 27, 2014 at 11:44 AM UTC
• Updated on: July 27, 2014 at 11:44 AM UTC
• Description:

  Attempts to use Contact Map information while compressing proteins in initial stages.

Code

--------------------------------------------------------------------
--           Code for Contact Mapper 
--------------------------------------------------------------------
-- by KarenCH
-- TODO:
--   Consider alternative underlying algorithms, such as SA or multiple starts or GAB
--   Consider mods to wigglers, such as having AllCMBands drop CI significantly during
--     squeeze.
--  What about a rebuild hooked up to one of the banders (disable - rebuild - enable - wiggle)
--  Try more actions such as "put one band" + rebuild at low ci; use ZLB-wiggler. "Til contact is 
--     made" banded wiggler. Put Comp-band modifications in.
--
-- If contact map info becomes "rich" again, with varying "heat" values, uncomment
-- the threshold slider in the dialog-handling function
--
--------------------------------------------------------------------

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

-- simple bander options
COREPULLING_CI = 0.70
BAND_STRENGTH_DEFAULT = 1.0   -- not really meaningful; wigglers strength-modify
BAND_STRENGTH_REGION = 1.0
BAND_MINCOMPRESS = 0.75
BAND_MAXEXPAND = 1.25
BANDEDWIGGLE_TRYCT = 5
BANDEDWIGGLE_TARGET = 500.0
BANDEDWIGGLE_SMALL_TARGET = 50.0
BIS_LENGTH = 5.0       -- max length of a BIS

-- values for controlling in-step wiggle timing and gain
QUICKWIGGLE_ITERS = 2
QUICKWIGGLE_GAIN  = 0.50    -- for iterative wiggles in quickwiggle time
PROB_QSTAB_EXTRA = 0.50

-- values for controlling behavior of 'fuse' (after an acceptable neighbor is found)
USE_SLOWWIGGLE = false   -- if true, then add touch-up fuse
SLOWWIGGLE_ITERS = 2
SLOWWIGGLE_GAIN    = 0.10
PROB_SLOWWIGGLE_CUTFUSE = 0.70
PROB_SLOWWIGGLE_LOCAL = 0.20
PROB_SLOWWIGGLE_BLUEFUSE = 0.10
SLOWWIGGLE_ONLYGLOBAL_THRESHOLD = 100.0

-- values for dealing with mutation
FORBIDMUTATION = true        -- I'm not playing with mutation at this time
MUTATE_NOT_SHAKE = false  -- this is for choosing between shake and mutate
NOMUTATELIST = {0}

-- values associated with idealize
IDEALIZE_USE_CUTS = true
IDEALIZE_FORBID_UNHEALED_CUTS = true
IDEALIZE_MAXTRIES = 5  -- a "try" fails if cuts not healed. increase ZLBstrength, try again.
IDEALIZE_ZLB_STRENGTH = 2.0
IDEALIZE_TUBE_RADIUS = 4.0
IDEALIZE_MAX_SEGLEN = 5

-- values associated with rebuilding
REBUILD_FAILSBEFOREQUIT = 75
REBUILD_TRYCOUNT = 5
REBUILD_ADDIDEALIZE = true
REBUILD_USECUTS = true
REBUILD_MAX_SEGLEN = 9
REBUILD_MAX_EARLY_SCOREDROP = 3000.0
REBUILD_MAX_MEDIUM_SCOREDROP = 500.0
REBUILD_MAX_LATE_SCOREDROP = 100.0

-- Clash Importance
SCALE_MAXCI = 1.0

-- values for contact map reading 
CONTACTMAP_THRESHOLD = 0.96
CONTACTMAP_ALLOW_TIPBANDS = true
CONTACTMAP_ADD_RANDOMBANDS = true
RANDOMBAND_COUNT = 3
CONTACTMAP_RESTRICTIONS = true     -- only improvements that don't worsen contact map

-------------------------------------------------------------------------------------------------------------
------------ VALUES USERS SHOULD NOT BE TOUCHING ---------------------------------
-------------------------------------------------------------------------------------------------------------

-------------- CONSTANTS THAT PROBABLY SHOULD NEVER CHANGE --------------
REBUILD_INNER_TUBE_RADIUS  = 12.0  -- 9.0 has been seen in some scripts
REBUILD_OUTER_TUBE_RADIUS = 18.0 
HELIX_DISTANCE_SKIP4 = 6.3  -- for skips of 4
HELIX_DISTANCE_SKIP5 = 8.6  -- for skips of 5

-- QUICKSAVE slots used herein  (QS_Run equivalent is taken by QS_TopLevelTemp)
QS_Start = 1
QS_Best = 3
QS_ShortUseTemp1 = 94   -- used by low-level subroutines that stash state for short time
QS_ShortUseTemp2 = 95 -- used by low-level subroutines that stash state for short time
QS_NeighborTemp1 = 96 -- used by "Neighbor Generation" subroutines (mid-level managers)
QS_NeighborTemp2 = 97 -- used by "Neighbor Generation" subroutines (mid-level managers)
QS_NeighborTemp3 = 98 -- used by "Neighbor Generation" subroutines (mid-level managers)
QS_TopLevelTemp = 99 -- used by top-level manager subroutines that call low-level managers

--------------- THINGS THE PROGRAM LEARNS AND REMEMBERS --------------
PuzzleAllowsMutate = false
PuzzleHasLockedSegs = false
PuzzleHasContactMap = false
IsDesignerPuzzle = false
EndCalled = false
InitialScore = 0.0
CurrentBestScore = 0.0
StartTime = 0
HasContactMapSegList = {}
ContactMapScore = 0.0
helixStarts = {}
helixEnds = {}
sheetStarts = {}
sheetEnds = {}
loopStarts = {}
loopEnds = {}
NonLockedSegList = {}
InitialClashImportance = behavior.GetClashImportance()
segCt = structure.GetCount()

---------------- STATISTICS-TRACKING ---------------------
algName = "none"
algNameList = {}
algStats = {
    name = "none",
    posScoreDelta = 0.0, 
    scoreDelta = 0.0, 
    runTime = 0, 
    callCount = 0,
    failCount = 0,
 }
TotalShakeTime = 0

------------------------------------------------------------------------------
--
---- CONTACT-MAPPER-SPECIFIC STUFF
--
------------------------------------------------------------------------------
OperationChosen = "none"

COUNT_OPERATIONS = 5000

----------------------------------------------------------------
-- HELPER FUNCTIONS
----------------------------------------------------------------
function getCMSegsWithinSphere( idxList, centerIdx, radius, includeCenter, heatThreshold )
    if includeCenter then
        idxList[ #idxList + 1 ] = centerIdx
    end
    for i=1, segCt do
        if i ~= centerIdx and structure.GetDistance( centerIdx, i ) <= radius and
            contactmap.GetHeat( centerIdx, i ) >= heatThreshold
        then
            idxList[ #idxList + 1 ] = i
        end
    end
end

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

----------------------------------------------------------------
-- BANDERS
----------------------------------------------------------------
function PutAllContactMapBands( heatThreshold, strength, doComp, doTip1, doTip2 )
    print( "AllCMBands" )
    changeSucceeded = false
    local baseStrength = strength
    for i = 1, segCt-2 do
        for j = i+2, segCt do
            if contactmap.GetHeat(i, j) >= heatThreshold then
                local goalLength = nil
                if doComp then
                    goalLength = structure.GetDistance( i, j )
                    if contactmap.IsContact( i, j ) then 
                        goalLength = 0.9 * goalLength
                        strength = 0.5 * baseStrength
                    else
                        goalLength = 0.7 * goalLength
                        strength = 1.5 * baseStrength
                    end
                else
                    if contactmap.IsContact( i, j ) then goalLength = structure.GetDistance( i, j ) end
                end
                changeSucceeded = BandBetweenSegsWithParameters( i, j, strength, goalLength, doTip1, doTip2 ) 
                                              or changeSucceeded
            end
        end
    end
    return changeSucceeded
end

function PutSomeContactMapBands( heatThreshold, strength, ctBands, doTip1, doTip2 )
    if ctBands > 1 then
        print( "SomeCMBands with ct="..ctBands )
    else
        print( "SingleCMBand" )
    end
    changeSucceeded = false
    local hotList = {}
    for i = 1, segCt-2 do
        for j = i+2, segCt do
            local heat = contactmap.GetHeat(i, j)
            if heat >= heatThreshold and not contactmap.IsContact( i , j ) then
                hotList[ #hotList + 1] = { i, j, heat }
            end
        end
    end

    randomizeIndexList( hotList )
    for i=1, math.min( ctBands, #hotList ) do
        --local which = random( #hotList )
        local ch = BandBetweenSegsWithParameters( hotList[i][1], hotList[i][2], strength, nil, doTip1, doTip2 ) 
        changeSucceeded = ch or changeSucceeded
    end
    return changeSucceeded
end

function PutAllContactMapBandsToSeg( idx, heatThreshold, strength, doComp, doTip1, doTip2 )
    print( "AllCMBandsToSeg ".. idx )
    changeSucceeded = false
    local baseStrength = strength
    if doComp then
        -- COMP-style bands
        for i=1, segCt do
            if i < idx - 1 or i > idx + 1 then
                local heat = contactmap.GetHeat( i, idx )
                if heat >= heatThreshold then
                    local goalLength = structure.GetDistance( i, idx )
                    if contactmap.IsContact( i, idx ) then
                        strength = 0.5 * baseStrength
                        goalLength = 0.9 * goalLength
                    else
                        strength = 1.5 * baseStrength
                        goalLength = 0.7 * goalLength
                    end
                    local changeSuc = BandBetweenSegsWithParameters( 
                                                  idx, i, strength, goalLength, 
                                                  doTip1, doTip2 ) 
                    changeSucceeded = changeSuc or changeSucceeded
                end
            end
        end
    else
        -- only bands to non-contacting segs
        local hotList = {}
        for i = 1, segCt do
            if i < idx - 1 or i > idx + 1 then
                local heat = contactmap.GetHeat( i, idx )
                if heat >= heatThreshold and not contactmap.IsContact( i , idx ) then
                    hotList[ #hotList + 1] = i
                end
            end
        end

        for i=1, #hotList do
            if hotList[ i ] ~= idx then  -- not supposed to happen!
                local changeSuc = BandBetweenSegsWithParameters( 
                                                  idx, hotList[i], strength * contactmap.GetHeat( i, idx ), 
                                                  nil, doTip1, doTip2 ) 
                changeSucceeded = changeSuc or changeSucceeded
            end
        end
    end
    return changeSucceeded
end

function PutContactMapBandsBetweenRegions( 
                                       heatThreshold, 
                                       startIdx1, endIdx1, startIdx2, endIdx2, 
                                       strength, doComp, doTip1, doTip2 )
    print( "CMTwoRegions: ("..startIdx1..", "..endIdx1..") to ("..startIdx2..", "..endIdx2..")" )
    changeSucceeded = false
    local baseStrength = strength
    for i = startIdx1, endIdx1 do
        for j = startIdx2, endIdx2 do
            if contactmap.GetHeat( i, j ) >= heatThreshold then
                local goalLength = nil
                if doComp then
                    goalLength = structure.GetDistance( i, j )
                    if contactmap.IsContact( i, j ) then 
                        strength = 0.5 * baseStrength
                        goalLength = 0.9 * goalLength 
                    else
                        strength = 1.5 * baseStrength
                        goalLength = 0.7 * goalLength 
                    end
                end
                changeSucceeded = BandBetweenSegsWithParameters( i, j, strength, goalLength, doTip1, doTip2 ) 
                                              or changeSucceeded
            end
        end
    end
    return changeSucceeded
end

function PutPushPullContactMapBands( badSeg, bandsWanted, center, strength, heatThreshold, wantPushNotPull )
    print( "PushPullCM: "..badSeg.." center="..center )
    changeSucceeded = false

    local okSegList = {}
    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[#okSegList + 1] = i
        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], strength, SR ) then 
            changeSucceeded = true 
            bandsAdded = bandsAdded + 1
        end
        if bandsAdded >= bandsWanted then break end
    end --  FOR(i)
    return changeSucceeded
end

function PutVoidKillerContactMapBands( idx, strength, heatThreshold )
    print( "VoidKillerCM at "..idx )
    changeSucceeded = false
    local t={} -- store possible segments

    for b=1,segCt do --test all segments
        if math.abs(idx-b) >= 15 and contactmap.GetHeat( idx, b ) >= heatThreshold then
            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, strength, d ) then changeSucceeded = true end
        end
    end

    return changeSucceeded
end

function PutCoreCompressBands(coreList, skipCt, startOffset, strength, compression, pullPhilics, doTips)
    if doTips == nil then doTips = false end
    if pullPhilics == nil then pullPhilics = false end
    print("CoreCompress: skipCt="..skipCt..", "..#coreList.." cores, and ".. TrimNum(compression) .. " c/e")
    local str = "   connecting to "
    if doTips then str=str .. "tips" else str = str .. "backbone" end
    if pullPhilics then str=str .. ", pulling philics away" else str = str .. ", philics and phobics treated same" end
    print( str )

    local maxToBand = #coreList
    local bandList = {}
    local added = 0

    -- 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
                local goalLength= compression * structure.GetDistance( coreList[ i ], coreList[ j ] )
                if BandBetweenSegsWithParameters( coreList[ i ], coreList[ j ], strength, goalLength ) then
                    added = added + 1
                end
            end
        end
    end

    -- now put bands between the "rest" of the protein and our chosen ones
    for i=startOffset, segCt, skipCt do
        local skipThisI = false
        for j=1, maxToBand do 
            if coreList [ j ] == i then 
                skipThisI = true
                break
            end
        end

        if not skipThisI then 
            for j=1, maxToBand do
                if i < coreList[ j ] - 1 or i > coreList[ j ] + 1 then
                    goalLength= compression * structure.GetDistance( i, coreList[ j ] )
                    if pullPhilics and not structure.IsHydrophobic( i ) then
                        goalLength= goalLength + 1.0
                    end
                    if BandBetweenSegsWithParameters( i, coreList[ j ], strength, goalLength, false, doTips ) then
                        added = added + 1
                    end
                end
            end
        end
    end

    return added > 0
end

function PutGentleCompBands(bandsWanted, minLength, maxLength, minSeparation, strength, compression )
    bandsAdded = 0
    failCount = 0
    repeat
        s = randomMovableSeg()
        e = randomMovableFarSeg( s, minSeparation, minLength, maxLength )
        if e == 0 then 
            failCount = failCount + 1
        else
            -- yes, this would try to add the same segment twice, but that will fail simply...
            changeSucceeded = BandBetweenSegsWithParameters( s,e, strength, dist * compression )
            if changeSucceeded then 
                bandsAdded = bandsAdded + 1
            else
                failCount = failCount + 1
            end
        end
    until bandsAdded == bandsWanted or failCount >= 100
    return true
end


----------------------------------------------------------------
-- ACTIONS
----------------------------------------------------------------
function PickSingleContactMapBand( )
    local iters = QUICKWIGGLE_ITERS
    changeSucceeded = false
    local heatThreshold = CONTACTMAP_THRESHOLD
    local strength = BAND_STRENGTH_DEFAULT
    local doTip1 = false
    local doTip2 = false
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTip1 = random() < 0.50
        doTip2 = random() < 0.50
    end
    local ctRandBands = 0
    if CONTACTMAP_ADD_RANDOMBANDS then
        if random( ) < 0.50 then ctRandBands = RANDOMBAND_COUNT end
    end

    changeSucceeded = PutSomeContactMapBands( heatThreshold, strength, 1, doTip1, doTip2 )
    if changeSucceeded then
        if ctRandBands > 0 then
            PutGentleCompBands( ctRandBands, 14.0, 30.0, 5, strength , 0.7 )
        end
        SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI ) 
    end
    return changeSucceeded
end

function PickSomeContactMapBands( )
    local iters = QUICKWIGGLE_ITERS
    changeSucceeded = false
    local heatThreshold = CONTACTMAP_THRESHOLD
    local strength = BAND_STRENGTH_DEFAULT
    local doTip1 = false
    local doTip2 = false
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTip1 = random() < 0.50
        doTip2 = random() < 0.50
    end
    local ctRandBands = 0
    if CONTACTMAP_ADD_RANDOMBANDS then
        if random( ) < 0.50 then ctRandBands = RANDOMBAND_COUNT end
    end
    local ctWanted = random( 5, 10 )

    changeSucceeded = PutSomeContactMapBands( heatThreshold, strength, ctWanted, doTip1, doTip2 )
    if changeSucceeded then
        if ctRandBands > 0 then
            PutGentleCompBands( ctRandBands, 14.0, 30.0, 5, strength , 0.7 )
        end
        SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI ) 
    end
    return changeSucceeded
end

function PickContactMapBandsForSeg( idx )
    if idx == nil then idx = randomContactableSeg( ) end
    changeSucceeded = false
    local iters = QUICKWIGGLE_ITERS
    local heatThreshold = CONTACTMAP_THRESHOLD
    local strength = BAND_STRENGTH_DEFAULT
    local doTip1 = false
    local doTip2 = false
    local doComp = random( ) < 0.50
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTip1 = random() < 0.50
        doTip2 = random() < 0.50
    end
    local ctRandBands = 0
    if CONTACTMAP_ADD_RANDOMBANDS then
        if random( ) < 0.50 then ctRandBands = RANDOMBAND_COUNT end
    end

    changeSucceeded = PutAllContactMapBandsToSeg( idx, heatThreshold, strength, doComp, doTip1, doTip2 )
    if changeSucceeded then
        if ctRandBands > 0 then
            PutGentleCompBands( ctRandBands, 14.0, 30.0, 5, strength , 0.7 )
        end
        if doComp then
            SimpleWiggleBeforeAndAfter( iters, 1, segCt, 1.0, false )
       else
            SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI ) 
       end
    end
    return changeSucceeded
end

function PickCMTwoRegionBander( idx1 )
    if idx1 == nil then idx1 = randomContactableSeg( ) end
    changeSucceeded = false
    local iters = QUICKWIGGLE_ITERS
    local strength = BAND_STRENGTH_DEFAULT
    local heatThreshold = CONTACTMAP_THRESHOLD
    local doTip1 = false
    local doTip2 = false
    local doComp = random( ) < 0.50
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTip1 = random() < 0.50
        doTip2 = random() < 0.50
    end
    local ctRandBands = 0
    if CONTACTMAP_ADD_RANDOMBANDS then
        if random( ) < 0.50 then ctRandBands = RANDOMBAND_COUNT end
    end

    local startIdx1, endIdx1 = GetRegionStartAndEnd( idx1 ) 
    if startIdx1 == 1 and endIdx1 == segCt then return false end -- we just can't deal with this
    local idx2 = startIdx1
    local tryCount = 0
    while idx2 >= startIdx1 and idx2 <= endIdx1 and tryCount < 50 do
        for i=startIdx1, endIdx1 do
            idx2 = randomContactSeg( i, CONTACTMAP_THRESHOLD )
            if idx2 ~= 0 and (idx2 < startIdx1 or idx2 > endIdx1) then break end -- got one
        end
        tryCount = tryCount + 1
    end
    if idx2 == 0 then return false end -- we have failed
    if tryCount == 50 then return false end
    local startIdx2, endIdx2 = GetRegionStartAndEnd( idx2 )
    changeSucceeded = PutContactMapBandsBetweenRegions( heatThreshold, 
                                                                      startIdx1, endIdx1, startIdx2, endIdx2, 
                                                                      strength, doComp, 
                                                                      doTip1, doTip2 )

    if changeSucceeded then
        if ctRandBands > 0 then
            PutGentleCompBands( ctRandBands, 14.0, 30.0, 5, strength , 0.7 )
        end
        if doComp then
            SimpleWiggleBeforeAndAfter( iters, 1, segCt, 1.0, false )
       else
            SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI ) 
       end
    end
    return changeSucceeded
end

function PickContactMapper( )
    local iters = QUICKWIGGLE_ITERS
    changeSucceeded = false
    local threshold = CONTACTMAP_THRESHOLD
    local strength = BAND_STRENGTH_DEFAULT
    local doTip1 = false
    local doTip2 = false
    local doComp = random( ) < 0.50
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTip1 = random() < 0.50
        doTip2 = random() < 0.50
    end
    local ctRandBands = 0
    if CONTACTMAP_ADD_RANDOMBANDS then
        if random( ) < 0.50 then ctRandBands = RANDOMBAND_COUNT end
    end

    changeSucceeded = PutAllContactMapBands( threshold, strength, doComp, doTip1, doTip2 )
    if changeSucceeded then
        if ctRandBands > 0 then
            PutGentleCompBands( ctRandBands, 14.0, 30.0, 5, strength , 0.7 )
        end
        if doComp then
            SimpleWiggleBeforeAndAfter( iters, 1, segCt, 1.0, false )
       else
            SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI ) 
       end
    end
    return changeSucceeded
end

function PickPushPullContactMapAction( badSeg )
    if badSeg == nil then badSeg = randomContactableSeg( ) end
    local segList = {}
    local iters = QUICKWIGGLE_ITERS
    local wantPushNotPull = random() < 0.50
    local ctBands = 5
    changeSucceeded = false
    local strength = BAND_STRENGTH_DEFAULT

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

    changeSucceeded = PutPushPullContactMapBands( badSeg, ctBands, centerSeg, strength, 
                                                                                  CONTACTMAP_THRESHOLD, wantPushNotPull )
    if changeSucceeded then 
        SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET ) 
    end
    
    DelBands()
    return changeSucceeded
end

function PickVoidKillerContactMapAction( idx )
    if idx == nil then idx = randomContactableSeg( ) end
    local iters = QUICKWIGGLE_ITERS
    local strength = BAND_STRENGTH_DEFAULT

    DelBands()
    changeSucceeded = PutVoidKillerContactMapBands( idx, strength, CONTACTMAP_THRESHOLD )
    if changeSucceeded then 
        SimpleBandedWiggleWrapper( iters, BANDEDWIGGLE_TARGET ) 
    end

    DelBands()
    return changeSucceeded
end

function PickCoreCompression( compression )
    local thisCore = {}
    local maxCenters = randomDice( 3, 1, 3 ) -- range 3 to 9, prefers 5-7
    local strength = BAND_STRENGTH_DEFAULT
    local skipCt = random( 3, 13 )
    local startOffset = random( 1, skipCt )
    local pullPhilics = random( ) < 0.5
    local doTips = false
    if CONTACTMAP_ALLOW_TIPBANDS then
        doTips = random() < 0.50
    end
 
    getCentralHydrophobics( thisCore, maxCenters)
     return PutCoreCompressBands( thisCore, skipCt, startOffset, strength, compression, pullPhilics, doTips )
end

function PickCoreSpreader( )
    local iters = QUICKWIGGLE_ITERS
    local maxBWigTries = 5
    changeSucceeded = false

    compression = random( 1.05, BAND_MAXEXPAND )
    local coreSegs = { }
    getCentralHydrophobics( coreSegs, math.floor(segCt / 8) )   -- for a 100-seg protein, this is about 13
 
    changeSucceeded = PickCoreCompression( compression )
    if not changeSucceeded then return false end

    BandedWiggleUntilEnoughChange( iters, BANDEDWIGGLE_TARGET, COREPULLING_CI )
    DelBands()
    return true
end

function PickStirThenContactMap( )
    PickCoreSpreader( )
    return PickContactMapper( )
end

----------------------------------------------------------------
-- THE CLEANUP WIGGLER
----------------------------------------------------------------
function DoCleanupWiggle( iters, minGainForIteration, prevScore )
    SetCI( 1.0 )
    local currScore = getScore( )

    recentbest.Save( )
    ShakeOrMutate( iters/2 )
    recentbest.Restore( ) -- just in case lowCI dropped the score instead of increasing it

    IterativeWiggleAll( iters, minGainForIteration, true, true)
    recentbest.Restore( ) -- just in case lowCI dropped the score instead of increasing it

    if USE_SLOWWIGGLE and prevScore - currScore < SLOWWIGGLE_ONLYGLOBAL_THRESHOLD  
    then
        local r = random( ) 
        if r < PROB_SLOWWIGGLE_LOCAL then
            print( "IterativeWiggleLocalByChunk" )
            IterativeWiggleLocalByChunk( iters, 1, segCt, minGainForIteration, random(3,6), random()<0.50 )
        elseif r < PROB_SLOWWIGGLE_LOCAL + PROB_SLOWWIGGLE_BLUEFUSE then
            print( "BlueFuse" )
            BlueFuse( )
        else   -- PROB_SLOWWIGGLE_CUTFUSE
            print( "IterativeCutAndWiggle" )
            IterativeCutAndWiggle( iters, minGainForIteration, random(3,11) )
        end
    end
    recentbest.Restore( ) -- just in case a fusish wiggle dropped the score instead of increasing it
end

----------------------------------------------------------------
-- CONTROLLERS and DIALOG BOXES
----------------------------------------------------------------
function PerformOperation(  )
    local minGainForIteration = SLOWWIGGLE_GAIN
    local iters = SLOWWIGGLE_ITERS

    local CtOperations = COUNT_OPERATIONS
    print( "Initial Score = "..TrimNum(getScore( )).." ("..GetContactScore( )..")  Performing " .. CtOperations .. " steps" )
    local OperationList = { 
                 "SimpleContactMapper", 
                 "StirThenContactMap", 
                 "SingleContactMapBand",
                 "FewContactMapBands",
                 "ContactMapBandsToSeg", 
                 "ContactMapBandTwoRegions",
                 "PushPullContactMap", 
                 "VoidKillerContactMap",
                 }

    algName = OperationChosen
    save.Quicksave( QS_Best )
    for i=1, CtOperations do
        save.Quickload( QS_Best )
        local startScore = getScore( )
        local cmScore = GetContactScore( )
        local startTime = os.time( )
        recentbest.Save( )

        if OperationChosen == "MixAndMatch" then
            algName = ChooseAlgorithm( OperationList )
        end

        print( "###Step "..i.." Score="..TrimNum(startScore).."  CmScore="..cmScore )
        local changeSucceeded = false
        if algName == "SimpleContactMapper" then
            changeSucceeded = PickContactMapper( )
        elseif algName == "StirThenContactMap" then
            changeSucceeded = PickStirThenContactMap( )
        elseif algName == "SingleContactMapBand" then
            changeSucceeded = PickSingleContactMapBand( )
        elseif algName == "FewContactMapBands" then
            changeSucceeded = PickSomeContactMapBands( )
        elseif algName == "ContactMapBandsToSeg" then
            changeSucceeded = PickContactMapBandsForSeg( )
        elseif algName == "ContactMapBandTwoRegions" then
            changeSucceeded = PickCMTwoRegionBander( )
        elseif algName == "PushPullContactMap" then
            changeSucceeded = PickPushPullContactMapAction( )
        elseif algName == "VoidKillerContactMap" then
            changeSucceeded = PickVoidKillerContactMapAction( )
        end

        -- Note: The operations can change algName by adding suffixes
        if not changeSucceeded then 
            print( "  Action Failed" )
            AddFailToStats()
        else
            local newScore = getScore( )
            local newCmScore = GetContactScore( )
            if USE_SLOWWIGGLE and newCmScore >= cmScore and newScore >= 0.995 * startScore then
                -- we're close... let's try again to get a gain...  was BlueFuse( )
                print( "CutAndWiggle because we're close" )
                IterativeCutAndWiggle( iters, minGainForIteration, random( 3,11 ) )
                didFuse = true
            end

            DelBands( )
            newScore = getScore( )
            newCmScore = GetContactScore( )
            local endTime = os.time( )
            print( "  Neighbor score: "..TrimNum( newScore ).." ("..newCmScore..")  time="..endTime - startTime )
            AddActionToStats( newScore - startScore, endTime - startTime )

            if newScore > startScore and (newCmScore >= cmScore or not CONTACTMAP_RESTRICTIONS) then
                print( ">>>> Action improved score to "..TrimNum(newScore).. "  CmScore="..newCmScore )

                -- yes, doing this can lower CM score, but we take it anyway (the broken points won't go far apart, at least)
                DoCleanupWiggle( iters, minGainForIteration, startScore )
                save.Quicksave( QS_Best )
            end
        end

        DelBands()
    end
end

function LetUserChooseParams( )
    dlg = dialog.CreateDialog("Contact Mapping Tricks")
    dlg.cmLabel1 = dialog.AddLabel( "----- Actions: Topmost checked box will be chosen-----" )
    dlg.cmMixAndMatch = dialog.AddCheckbox( "Mix and Match actions (all)", true )
    dlg.cmStirrerMap = dialog.AddCheckbox( "Stirred ContactMapper", false )
    dlg.cmSimpleMap = dialog.AddCheckbox( "Full ContactMapper", false )
    dlg.cmOneMapBand = dialog.AddCheckbox( "Single CM band", false )
    dlg.cmFewMapBands = dialog.AddCheckbox( "Few scattered CM bands", false )
    dlg.cmMapBandsToSeg = dialog.AddCheckbox( "All CM bands to a seg", false )
    dlg.cmMapBandsTwoRegions = dialog.AddCheckbox( "CM bands between two regions", false )
    dlg.cmPushPull = dialog.AddCheckbox( "PushPull with ContactMap", false )
    dlg.cmVoidKiller = dialog.AddCheckbox( "VoidKiller with ContactMap", false )
    dlg.cmLabel2 = dialog.AddLabel( "--------------------- Options -------------------" )
    dlg.cmOpCount = dialog.AddSlider( "Total Tries", COUNT_OPERATIONS, 100, 100000, 0 )
    dlg.cmMaxCI = dialog.AddSlider( "Max CI", SCALE_MAXCI, 0.00, 1.00, 2 )
    -- dlg.cmHeatThreshold = dialog.AddSlider( "HeatMinThreshold", CONTACTMAP_THRESHOLD, 0.0, 1.0, 2 )
    dlg.cmAllowTipBands = dialog.AddCheckbox( "Allow bands to tips", CONTACTMAP_ALLOW_TIPBANDS )
    dlg.cmAddRandBands = dialog.AddCheckbox( "Add random bands", CONTACTMAP_ADD_RANDOMBANDS )
    dlg.cmForceContactBettering = dialog.AddCheckbox( "Contact map tries to improve", CONTACTMAP_RESTRICTIONS )
    dlg.cmFuseSuccesses = dialog.AddCheckbox( "Follow-on Fuse", USE_SLOWWIGGLE )

    dlg.ok = dialog.AddButton("OK", 1)
    dlg.cancel = dialog.AddButton("Cancel", 0)
    if dialog.Show(dlg) == 0 then return false end

    if dlg.cmMixAndMatch.value then OperationChosen = "MixAndMatch"
    elseif dlg.cmStirrerMap.value then OperationChosen = "StirThenContactMap"
    elseif dlg.cmSimpleMap.value then OperationChosen = "SimpleContactMapper"
    elseif dlg.cmOneMapBand.value then OperationChosen = "SingleContactMapBand"
    elseif dlg.cmFewMapBands.value then OperationChosen = "FewContactMapBands"
    elseif dlg.cmMapBandsToSeg.value then OperationChosen = "ContactMapBandsToSeg"
    elseif dlg.cmMapBandsTwoRegions.value then OperationChosen = "ContactMapBandTwoRegions"
    elseif dlg.cmPushPull.value then OperationChosen = "PushPullContactMap"
    elseif dlg.cmVoidKiller.value then OperationChosen = "VoidKillerContactMap"
    else return false -- something went wrong...
    end

    COUNT_OPERATIONS = dlg.cmOpCount.value
    --CONTACTMAP_THRESHOLD = dlg.cmHeatThreshold.value
    CONTACTMAP_RESTRICTIONS = dlg.cmForceContactBettering.value
    USE_SLOWWIGGLE = dlg.cmFuseSuccesses.value
    CONTACTMAP_ALLOW_TIPBANDS = dlg.cmAllowTipBands.value
    CONTACTMAP_ADD_RANDOMBANDS = dlg.cmAddRandBands.value
    SCALE_MAXCI = dlg.cmMaxCI.value

    if CONTACTMAP_RESTRICTIONS then 
        print( "Contact map scores may not drop" ) 
    else 
        print( "Contact scores can drop" ) 
    end
    if USE_SLOWWIGGLE then 
        print( "Touch-up fusing will happen" ) 
    else 
        print( "No fuses" ) 
    end
    if CONTACTMAP_ALLOW_TIPBANDS then 
        print( "Bands can connect to sidechain tips" ) 
    else 
        print( "Bands only to backbone" ) 
    end
    if CONTACTMAP_ADD_RANDOMBANDS then 
        print( "Can add ".. RANDOMBAND_COUNT.. " random extra bands" ) 
    else 
        print( "No random extra bands" ) 
    end
    return true
end


----------------------------------------------------------------
-- 
----------------    BOILERPLATE    ---------------------
-- 
----------------------------------------------------------------

----------------------------------------------------------------
-- BASIC FUNCTIONALITY
----------------------------------------------------------------
function DebugPrint( str )
    if DEBUGRUN then print( str ) end
end

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

-- Not for "external" use - call getScore. This could change if customers want 
-- something besides current or recentbest.
function internalGetScore( 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 getScore( )
    return internalGetScore( false )
end
function getRBScore( )
    return internalGetScore( true )
end

function SaveBest( )
    local score = getScore( )
    if score > CurrentBestScore then
        save.Quicksave( QS_Best )
        CurrentBestScore = score
    end
end

function LoadBest( )
    save.Quickload( QS_Best )
end

function SetCI( ci )
    behavior.SetClashImportance( SCALE_MAXCI * ci )
end

function ExactSetCI( ci )
    behavior.SetClashImportance( ci )
end

----------------------- MATHY STUFF -----------------------
function seedRandom()
    seed=os.time( )/math.abs( getScore( ) )
    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 )

    -- throw away a couple of randoms
    math.random( )
    math.random( )
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 == 0 then return 0 end -- a random number between 0 and 0 is 0
            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 randomMovableFarSeg( segOrigin, minSep, minDist, maxDist )
    local segList = {}
    for e=1, segCt do
        if IsMovableSeg( e ) and (e <= segOrigin - minSep or e >= segOrigin + minSep ) then
            dist = structure.GetDistance( segOrigin, e )
            if dist >= minDist and dist <= maxDist then segList[ #segList + 1 ] = e end
        end
    end
    if #segList >= 1 then return segList[ random( #segList ) ]
    else return 0
    end
end

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

function randomContactableSeg( )
    return HasContactMapSegList[ random( #HasContactMapSegList ) ]
end

function randomContactSeg( segIn, heatThreshold )
    local segOpts = {}
    for i = 1, segCt do
        if i < segIn - 1 or i > segIn + 1 then
            if contactmap.GetHeat( segIn, i ) >= heatThreshold then
                segOpts[ #segOpts + 1] = i
            end
        end
    end
    if #segOpts == 0 then return 0 end
    return segOpts[ random( #segOpts ) ]
end

function RandomPrimeNotMoreThan( max )
    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}
    lastIdx = 1
    for i=1, #primes do
        if primes[ i ] > max then break end
        lastIdx = i
    end
    return primes[ math.random( lastIdx ) ]
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

function ConstructScoreMatrix( scmat )
    -- TODO: Get the list of "active" subscores from puzzle info; use that info here.
    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 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

---------- 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

function SortIdxListBySubscore( idxList, subscoreName )
    local subscoreList = {}
    for i=1, #idxList do
        subscoreList[ i ] =  current.GetSegmentEnergySubscore( idxList[ i ], subscoreName )
    end
    sortAllByScore( idxList, subscoreList )
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)
    for i=1, #idxListInternal do
        if structure.IsHydrophobic( idxListInternal[ i ] ) then
            idxList[ #idxList + 1 ] = idxListInternal[ i ]
            if #idxList == ctWanted then return end
        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

----------------------------------------------------------------
-- MULTI-TRY ACTION CONTROLLER
----------------------------------------------------------------
function MultitryActOnRegion( NumberOfStarts, NumberOfActionsPerTry, startIdx, endIdx, ActionFunction )
    print( "MultitryActOnRegion for "..NumberOfStarts.." starts with "..NumberOfActionsPerTry.." steps" )
    BestScore=current.GetScore( )
    LoadBest( )

    save.Quicksave( QS_NeighborTemp1 )
    save.Quicksave( QS_NeighborTemp2 )
    save.Quicksave( QS_NeighborTemp3 )

    for k=1, NumberOfStarts do
        save.Quickload( QS_NeighborTemp1 )   -- always start over from initial state
        save.Quicksave( QS_NeighborTemp2 )
        print( "Starting try " .. k.. ": Score="..getScore( ) )

        local initscore = getScore( )
        for j=1, NumberOfActionsPerTry do     -- use "best so far"
            save.Quickload( QS_NeighborTemp2 )
            local currScore = getScore( )

            ActionFunction( startIdx, endIdx )

            local endScore = getScore( )
            if endScore > currScore then
                print( ">>>> At ".. j ..": score="..endScore )
                save.Quicksave( QS_NeighborTemp2 )
                currScore = endScore
            end
            if endScore > BestScore then
                save.Quicksave( QS_NeighborTemp3 )
                BestScore = endScore
            end
            SaveBest( )
        end    -- FOR j<actions

        save.Quickload( QS_NeighborTemp2 )
        local endScore = getScore( )
        print( "Finishing try " .. k.. "  Score: "..endScore  )
    end       -- FOR k<starts

    LoadBest( )
end

----------------------------------------------------------------
-- ACTION STATISTICS
----------------------------------------------------------------
function WhereIsStringInArray( array, str )
    for i=1, #array do if str == array[ i ] then return i end end
    return 0
end

function AddActionToStats( scoreDelta, timeDelta )
    if WhereIsStringInArray( algNameList, algName ) == 0 then
        algNameList[ #algNameList + 1] = algName
    end

    if algStats[ algName ] == nil then 
        algStats[ algName ] = { 
                  posScoreDelta = 0.0, 
                  scoreDelta = 0.0, 
                  runTime = 0, 
                  callCount = 0,
                  gainCount = 0,
                  failCount = 0,
        }
    end

    if scoreDelta > 0.0 then
        algStats[ algName ].posScoreDelta = algStats[ algName ].posScoreDelta + scoreDelta
        algStats[ algName ].gainCount = algStats[ algName ].gainCount + 1
    end
    algStats[ algName ].scoreDelta = algStats[ algName ].scoreDelta + scoreDelta
    algStats[ algName ].runTime = algStats[ algName ].runTime + timeDelta
    algStats[ algName ].callCount = algStats[ algName ].callCount + 1
end

function AddFailToStats( )
    if WhereIsStringInArray( algNameList, algName ) == 0 then
        algNameList[ #algNameList + 1] = algName
    end

    if algStats[ algName ] == nil then 
        algStats[ algName ] = { 
                  posScoreDelta = 0.0, 
                  scoreDelta = 0.0, 
                  runTime = 0, 
                  callCount = 0,
                  gainCount = 0,
                  failCount = 1,
        }
        return
    end

    algStats[ algName ].failCount = algStats[ algName ].failCount + 1
end

function ChooseAlgorithm( algList )
    -- we dither runs and gains by 1 partly to avoid divide-by-zero, and 
    -- partly to encourage bringing in not-yet-used algorithms.
    local scaleValue = 0.0
    local funcValues = {}

    for i = 1, #algList do
        local runs = 0
        local gains = 0
        if algStats[ algList[ i ] ] == nil then 
            runs = 1
            gains = 1
        else
            -- failCount/2 is a bit artificial, but failed actions should count against an algorithm too.
            runs = 1 + algStats[ algList[ i ] ].callCount + math.floor( algStats[ algList[ i ] ].failCount / 2 )
            gains = 1 + algStats[ algList[ i ] ].gainCount
        end
        funcValues[ i ] = gains / runs    -- finding value of this item
        scaleValue = scaleValue + funcValues[ i ]   -- finding total size of bucket
    end

    local which = random( )
    local sum = 0.0
    for i = 1, #funcValues do
        sum = sum + ( funcValues[ i ] / scaleValue )
        if which < sum then 
            return algList[ i ] 
        end
    end
    algName = algList[ #algList ]
    return algName
end

function ReportStats( )
    print("Action stats:")
    for i = 1, #algNameList do
        local t = algNameList[ i ]
        if algStats[ t ] ~= nil and algStats[ t ].callCount > 0 then
            print( "  " .. t .. ": scoreDel=" .. TrimNum( algStats[ t ].scoreDelta ) ..
                     " posScore=" .. TrimNum( algStats[ t ].posScoreDelta ) ..
                     " runTime=" .. algStats[ t ].runTime ..
                     " runs=" .. algStats[ t ].callCount ..
                     " gains=" .. algStats[ t ].gainCount ..
                     " fails=" .. algStats[ t ].failCount )
        end
    end
end

function ResetStats( )
    for i = 1, #algNameList do
        local t = algNameList[ i ]
        if algStats[ t ] ~= nil then
            algStats[ t ].scoreDelta = 0
            algStats[ t ].posScoreDelta = 0
            algStats[ t ].runTime = 0
            algStats[ t ].callCount = 0
            algStats[ t ].gainCount = 0
            algStats[ t ].failCount = 0
        end
    end
end

----------------------------------------------------------------
-- INFORMATION GATHERING
----------------------------------------------------------------

---------------- SIDECHAIN FUNCTIONS -------------------
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

-- for branched aas, simply picks one (longer, one with donor/acceptor tip, or if no difference then either)
function GetAtomOfTip( aa )
    if aa == "a" then return 10
    elseif aa == "c" then return 11
    elseif aa == "d" then return 8
    elseif aa == "e" then return 9
    elseif aa == "f" then return 20
    elseif aa == "g" then return 0 -- glycine has no tip. just use a backbone atom
    elseif aa == "h" then return 17
    elseif aa == "i" then return 18
    elseif aa == "k" then return 9
    elseif aa == "l" then return 16
    elseif aa == "m" then return 17
    elseif aa == "n" then return 14
    elseif aa == "p" then return 13
    elseif aa == "q" then return 9
    elseif aa == "r" then return 22
    elseif aa == "s" then return 11
    elseif aa == "t" then return 6
    elseif aa == "v" then return 13
    elseif aa == "w" then return 24
    elseif aa == "y" then return 12
    else return 0
    end
end

function GetTipAtomOfSeg( idx )
    return GetAtomOfTip( structure.GetAminoAcid( idx ))
end

-------------- DESIGNER PUZZLE PROPERTIES  --------------
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

------------- LOCKED-SEGMENT-RELATED FUNCTIONS ------------
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

------------- CONTACTMAP FUNCTIONS ------------
function CheckForContactMap( )
    PuzzleHasContactMap = false
    local saveval = 0.0
    for i=1, segCt-1 do
        for j=i+1, segCt do
            val = contactmap.GetHeat( i, j )
            if saveval  ~= 0.0 and val ~= saveval then
                PuzzleHasContactMap = true
                ContactMapScore = GetContactScore( )
                return   -- all we wanted to know was whether scores exist
            end
            if saveval == 0.0 then saveval = val end
        end
    end
    return
end

function SegHasContactData( segIn, heatThreshold )
    for i = 1, segCt do
        if i < segIn - 1 or i > segIn + 1 then
            if contactmap.GetHeat( segIn, i ) >= heatThreshold then return true end
        end
    end
    return false
end

function InitializeContactMapSegList( heatThreshold )
    HasContactMapSegList = {}
    for i = 1, segCt do
         if SegHasContactData( i, heatThreshold ) then 
             HasContactMapSegList[ #HasContactMapSegList + 1 ] = i 
        end
    end
end

function GetContactScore( )
    if not PuzzleHasContactMap then return 0 end
    local sum = 0.0
    local segCt = structure.GetCount( )
    for i = 1,segCt-1 do
        for j = i + 1, segCt do
            if contactmap.IsContact( i, j ) then sum = sum + contactmap.GetHeat( i, j ) end
         end
    end
    return sum
end

function ActionDamagedContactMap( )
    if PuzzleHasContactMap then 
        local sc = GetContactScore( )
        if sc < ContactMapScore then return true end
        ContactMapScore = sc   -- now we have a "better" score, so keep that one
    end
    return false
end

------------- REGION-IDENTIFICATION FUNCTIONS ------------
function getHelices()
    inside = false
    for i=1, segCt do
        if ( IsMovableSeg( i ) and structure.GetSecondaryStructure(i) == "H" ) then
            if ( inside == false ) then
                inside = true
                helixStarts[ #helixStarts + 1] = i
            end
        elseif ( IsMovableSeg( i ) and 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 ( IsMovableSeg( i ) and structure.GetSecondaryStructure(i) == "E" ) then
            if ( inside == false ) then
                inside = true
                sheetStarts[ #sheetStarts + 1 ] = i
            end
        elseif (IsMovableSeg( i ) and  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 ( IsMovableSeg( i ) and structure.GetSecondaryStructure(i) == "L" ) then
            if ( inside == false ) then
                inside = true
                loopStarts[ #loopStarts + 1 ] = i
            end
        elseif ( IsMovableSeg( i ) and 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
    if not IsMovableSeg( idx ) then return idx,idx end

    local typeIdx = structure.GetSecondaryStructure( idx )
    if idx > 1 then
        s = idx - 1
        while s >= 1 do
            if structure.GetSecondaryStructure( s ) ~= typeIdx or not IsMovableSeg( s ) 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  or not IsMovableSeg( e ) 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)
    if not IsMovableSeg( idx ) then 
        return idx, idx 
    end
    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

----------------------------------------------------------------------
-- CUTTING FUNCTIONS 
----------------------------------------------------------------------
function DeleteCutAtIndex( 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 deletion didn't happen!
end

function AddRandomCuts( cutList, maxCuts )
    for i=1, maxCuts do
        local seg = randomMovableSeg( )
        local gotMatch = false
        for j=1, #cutList do 
            if cutList[ j ] == seg then gotMatch = true end 
        end
        if not gotMatch then
            structure.InsertCut( seg )
            cutList[ #cutList + 1 ] = seg
        end
    end
end

function CleanCutsInList( cutList )
    changeSucceeded = true
    for i=1, #cutList do
        changeSucceeded = DeleteCutAtIndex( cutList[ i ] ) and changeSucceeded
    end
    return changeSucceeded
end

----------------------------------------------------------------------
-- BARE-BONE REBUILD FUNCTION 
----------------------------------------------------------------------
function RebuildSelected( iters )
    band.DisableAll( )
    structure.RebuildSelected( iters )
    band.EnableAll( )
end
------------------------------------------------------------------------------
--
---- WIGGLER FUNCTIONS
--
------------------------------------------------------------------------------

------------------------------------------------------------------------------
-- GENERAL-PURPOSE SHAKERS AND WIGGLERS
------------------------------------------------------------------------------
function ShakeOrMutate( iters, idxList )
    local startTime = os.time( )
    if iters == nil then iters = 1 end
    if idxList == nil then
        if MUTATE_NOT_SHAKE and not FORBIDMUTATION then
            structure.MutateSidechainsAll( iters )
        else
            if not IsDesignerPuzzle then structure.ShakeSidechainsAll( iters ) end
        end
    else
        selection.DeselectAll( )
        for i=1, #idxList do selection.Select( i ) end
        if MUTATE_NOT_SHAKE and not FORBIDMUTATION then
            structure.MutateSidechainsSelected ( iters )
        else
            if not IsDesignerPuzzle then structure.ShakeSidechainsSelected ( iters ) end
        end
        selection.DeselectAll( )
    end
    TotalShakeTime = TotalShakeTime + ( os.time( ) - startTime )
end

function WiggleRange( iters, 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( iters, true, true )
    else
        structure.WiggleSelected( iters, 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 IterativeWiggleAll( iters, minGain, doBackbone, doSidechain )
    local lastScore = getScore( )
    local gain = minGain
    while gain >= minGain do
        structure.WiggleAll( iters, doBackbone, doSidechain )  -- do's ok for nil
        gain = getScore( ) - lastScore
        lastScore = getScore( )
    end
end

---------------------------------------------
-- SIMPLE STABILIZER WIGGLES
---------------------------------------------
function IterativeWiggleList( idxList, iters, 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( iters, doBackbone, doSidechain )  -- do's ok for nil
        gain = getScore( ) - lastScore
        lastScore = getScore( )
    end
    selection.DeselectAll( )
end

function IterativeWiggleShakeWiggle( idxList, iters )
    IterativeWiggleList( idxList, iters, QUICKWIGGLE_GAIN )
    sc = getScore( )
    ShakeOrMutate( 1, idxList )
    if math.abs( getScore( ) - sc) > 1 then    -- why bother if s/m didn't do anything?
        IterativeWiggleAll(iters, QUICKWIGGLE_GAIN, true, true )
    end
end

function qStabWiggle( iters, doExtraWork )
    if doExtraWork == nil then doExtraWork = false end
    SetCI( 1.0 )
    ShakeOrMutate( 1 )

    if doExtraWork then
        SetCI( 0.4 )
        IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
        SetCI( 1.0)
        ShakeOrMutate( 1 )
    end
    SetCI( 1.0)
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )   -- was 10*
end

-- the only user for now is the idealizer
function WiggleZlbNicely( iters, strength, skipLower, skipUpper )
    if skipLower == nil then skipLower = segCt+1 end
    if skipUpper == nil then skipUpper = 0 end
    ZLB(strength, skipLower, skipUpper)
    --structure.WiggleAll( iters, true, true)
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
end

function SimpleRegionWiggle( iters, idxLo, idxHi, pickGlobalWiggle, wiggleExpansion )
    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 )
    WiggleRange( iters, startIdx, endIdx, false, pickGlobalWiggle )
end

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

    SetCI( pullingCI )
    deltaPoints = 0.0
    ss=getScore()

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

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

-- returns the actual amount of change that occurred (absolute value) 
function BandedWiggleUntilTargetMet( iters, scoreDeltaTarget, maxTries, pullingCI )
    print( "BandedWiggleUntilTargetMet" )
    save.Quicksave( QS_ShortUseTemp )
    save.Quicksave( QS_ShortUseTemp2 )
    if pullingCI == nil then pullingCI = COREPULLING_CI end
    SetCI( pullingCI )

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

        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 
            save.Quicksave( QS_ShortUseTemp2 )
            break   -- 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...
    if getRBScore( ) > getScore( ) then 
        recentbest.Restore( )  -- sometimes score went up. farbeit from us to refuse an improvement
    end
    SetCI( 1.0 )

    DelBands( )
    qStabWiggle( iters, random() < PROB_QSTAB_EXTRA )  -- a bit of wiggle afterwards to help settle down
end

function SimpleWiggleBeforeAndAfter( iters, startIdx, endIdx, probQstab, forbidGlobalWiggle )
    print( "SimpleWiggleBeforeAndAfter" )
    if forbidGlobalWiggle == nil then forbidGlobalWiggle = false end
    local wantGlobalWiggle = false
    if not forbidGlobalWiggle then wantGlobalWiggle = random() < 0.50 end
    local pullingCI = COREPULLING_CI

    if random() < probQstab then
        qStabWiggle( iters, random() < PROB_QSTAB_EXTRA ) -- 2nd param is a "work harder" one
    else
        SetCI( pullingCI )
        SimpleRegionWiggle( iters, startIdx, endIdx, wantGlobalWiggle, random(4) )
        SetCI( 1.0 )
        IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
    end

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

    if random() < probQstab then
        qStabWiggle( iters, random() < PROB_QSTAB_EXTRA ) -- 2nd param is a "work harder" one
    else
        ShakeOrMutate( iters / 2 )
        wantGlobalWiggle = not wantGlobalWiggle   -- invert decision about global wiggles.
        if forbidGlobalWiggle then wantGlobalWiggle = false end
        SimpleRegionWiggle( iters, startIdx, endIdx, wantGlobalWiggle, random(4,8) )
        IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
    end
end

function nudge( score )
    SetCI( 0.3 )
    structure.WiggleSelected( 1 )
    SetCI( 1 )
    structure.WiggleSelected( 8 )
end
function WiggleAndShakeForRebuilder( initScore, maxFall )
    print( "WiggleAndShakeForRebuilder" )
     -- initScore is the score before our rebuild happened
    recentbest.Save( )

    scoreA = getScore( )
    if initScore - scoreA > maxFall then
        ShakeOrMutate( 2 )
    elseif initScore - scoreA > maxFall/2 then
        ShakeOrMutate( 1 )
    end

    scoreA = getScore( )
    structure.WiggleAll( 15 )
    recentbest.Restore( )

    scoreB = getScore( )
    if scoreB - scoreA < 10 and initScore - scoreB > 30 then
        -- The wiggle got stuck and didn't achieve anything. try a nudge to see if that helps
        nudge( )
        recentbest.Restore( )
        scoreB = getScore( )
    end

    if initScore - scoreB <= 100.0 then
        ShakeOrMutate( 1 )
        structure.WiggleSelected( 10 )
        recentbest.Restore( )
        nudge( )
        recentbest.Restore( )
    end   -- else not worth the effort
end

--------------------------------------------------
-- SLOW CLEANUP WIGGLERS
--------------------------------------------------
function BlueFuse( )
    recentbest.Save( )

    SetCI(.05)
    structure.ShakeSidechainsAll(1)
    SetCI(1)
    structure.WiggleAll(5)

    recentbest.Restore( )

    SetCI(.07)
    structure.ShakeSidechainsAll(1)
    SetCI(1)
    structure.WiggleAll(5)

    recentbest.Restore( )

    SetCI(.3)
    structure.WiggleAll(8)
    SetCI(1)
    structure.WiggleAll(15)

    recentbest.Restore( )
end

function IterativeWiggleLocalByChunk(iters, 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
            recentbest.Save( )
            startIdx = startSeg + currOffset + chunkSize*idx
            WiggleRange(iters, startIdx, startIdx + chunkSize, addFreeze, true)
            recentbest.Restore( )
            if addFreeze then freeze.UnfreezeAll() end
        end
        local currScore = getScore( )
        gain = currScore - lastScore
        lastScore = currScore
    end
end

function BasicCutAndWiggleFuse( baseIters, shift, offset )
    if shift == nil then shift = 3 end
    if offset == nil then offset = 1 else offset = 1 + (offset % shift) end

    for i = offset, segCt-1, shift do
        structure.InsertCut(i)
    end
    recentbest.Save( )  -- this save/restore is entirely within cut region, so no cut problems, probably
    WiggleWiggleSimpleFuse( 0.30, 1.0, baseIters, 2*baseIters )
    recentbest.Restore( )
    for i = 1, segCt do -- take the opportunity to kill any "lying around" cuts
        DeleteCutAtIndex( i )
    end

    recentbest.Save( )
    WiggleWiggleSimpleFuse( 0.30, 1.0, baseIters, 4*baseIters )
    recentbest.Restore( )

    SetCI( 1.0 )
end

function IterativeCutAndWiggle( baseIters, minGain, shift )
    save.Quicksave( QS_ShortUseTemp1 )

    local score = getScore( )
    local gain = minGain
    repeat
        save.Quickload( QS_ShortUseTemp1 )
        BasicCutAndWiggleFuse( baseIters, shift, random( shift ) )
        local newscore = getScore( )
        if newscore >= score then 
            gain = newscore - score
            score = newscore
            save.Quicksave( QS_ShortUseTemp1 )
        else
            gain = 0.0
        end
    until gain < minGain
    save.Quickload( QS_ShortUseTemp1 )
end

--------------------------------------------------
-- BASIC FUSE WIGGLERS
--------------------------------------------------
function WiggleShakeSimpleFuse( ci, iters )
    SetCI( ci )
    structure.WiggleAll( iters )  -- instead of IterativeWiggleAll...
    ShakeOrMutate( iters/2 )
end
function WiggleWiggleSimpleFuse( ci1, ci2, shortiter, longiter )
    SetCI( ci1 )
    structure.WiggleAll( shortiter )
    SetCI( ci2 )
    structure.WiggleAll( longiter )  -- "long" wiggle
end
function ShakeWiggleFuse( ci1, ci2, iters )  -- aka Fuze1
    SetCI( ci1 )
    ShakeOrMutate( iters/2 )
    SetCI( ci2 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
end
function WiggleWiggleWiggleFuse( ci1, ci2, iters ) -- aka Fuze2
    SetCI( ci1 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
    SetCI( 1.00 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
    SetCI( ci2 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
end
function WiggleShakeWiggleFuse( iters )  -- aka FuzeEnd
    SetCI( 1.00 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
    ShakeOrMutate( 1 )
    IterativeWiggleAll( iters, QUICKWIGGLE_GAIN, true, true )
end

----------------------------------------------------------------------
--
---- BANDER FUNCTIONS
--
----------------------------------------------------------------------
function DelBands()
    band.DeleteAll()
end

function BandBetweenSegsWithParameters( seg1, seg2, strength, goalLength, fromTip1, fromTip2 )
    -- caller can set goalLength without setting strength by providing strength = 0.0
    if not ( IsMovableSeg( seg1) or IsMovableSeg( seg2) ) then return false end
    if fromTip1 == nil then fromTip1 = false end
    if fromTip2 == nil then fromTip2 = false end
    if fromTip1 then tip1 = GetTipAtomOfSeg(seg1) else tip1 = 0 end
    if fromTip2 then tip2 = GetTipAtomOfSeg(seg2) else tip2 = 0 end

    local bIdx = band.AddBetweenSegments( seg1, seg2, tip1, tip2 )
    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

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

----------------------------------------------------------------------
--
---- IDEALIZERS AND REBUILDERS
--
----------------------------------------------------------------------

----------------------------------------------------------------------
-- IDEALIZERS
----------------------------------------------------------------------
function IdealizeInRange( startSeg, endSeg, iters, useCuts )
    print( "IdealizeInRange for "..startSeg.." to "..endSeg )
    if useCuts == nil then useCuts = true end
    local idxList = {}
    changeFailed = false
    if startSeg > endSeg then startSeg, endSeg = endSeg, startSeg end
    if startSeg < 1 then startSeg = 1 end
    if endSeg > segCt then endSeg = segCt end
    local tries = 0

    save.Quicksave( QS_ShortUseTemp1 )

    local delta = 0.0
    repeat
        save.Quickload( QS_ShortUseTemp1 )
        changeFailed = false
        if IDEALIZE_MAXTRIES > 1 and tries == IDEALIZE_MAXTRIES - 1 then
            delta = 0.0
            useCuts = false   -- overrule the user's request for use of cuts (should always succeed)
        end

        if useCuts then
            if startSeg - 2 > 1 then 
                DebugPrint( "  putting cut at ".. startSeg - 2 )
                structure.InsertCut( startSeg - 2 ) 
            end
            if endSeg + 1 < segCt then 
                DebugPrint( "  putting cut at ".. endSeg + 1 )
                structure.InsertCut( endSeg + 1 ) 
            end

            if startSeg - 3 >= 1 then
                BandBetweenSegsWithParameters( startSeg-3, startSeg, 1.0 + delta, 
                                                                      structure.GetDistance( startSeg-3, startSeg ) )
            end

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

        selection.DeselectAll()
        getSegsWithinRangeOfTube( idxList, startSeg, endSeg, IDEALIZE_TUBE_RADIUS )
        for i=1, #idxList do selection.Select( idxList[i] ) end
        structure.IdealizeSelected()
        selection.DeselectAll()

        if useCuts then
            -- get some wiggling happening
            print( "WiggleZlbNicely" )
            WiggleZlbNicely( iters, IDEALIZE_ZLB_STRENGTH + delta, 
                                      math.max( 1, startSeg-2 ), math.min( segCt, endSeg+2 ) )

            -- now put the protein back in its proper shape
            if startSeg - 2 > 1 then
                if not DeleteCutAtIndex(startSeg - 2 ) and IDEALIZE_FORBID_UNHEALED_CUTS then 
                    changeFailed = true 
                end
            end
            if changeFailed == false and endSeg + 1 < segCt then 
                if not DeleteCutAtIndex(endSeg + 1 ) and IDEALIZE_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
            print( "WiggleRange" )
            WiggleRange( iters, startSeg - 2, endSeg + 1, false, random() < 0.50 )
            ShakeOrMutate( iters/2 )
        end

        if not changeFailed then 
            DelBands()
            return true 
        end
        tries = tries + 1
        delta = delta + 0.25
    until tries >= IDEALIZE_MAXTRIES
    DelBands()
    return false
end

----------------------------------------------------------------------
-- REBUILDERS
-- TRY: freeze before and after start-end or cut before and after start-end
-- ALSO: try a rebuild that changes Secondary structure to loop before
-- run, restores it after run
----------------------------------------------------------------------
function RebuildInRange( startSeg, endSeg, iters, maxTries, maxFall )
    if maxFall == nil then maxFall = REBUILD_MAX_MEDIUM_SCOREDROP end
    print( "RebuildInRange from "..startSeg.." to "..endSeg )
    changeSucceeded = false
    local oldScoreMat = {}
    local locss = {}
    local killedSS = false
    local startScore = getScore( )

    ConstructScoreMatrix( oldScoreMat )

    -- probably PickRebuild should decide whether to do secondary structure stuff
    for i=startSeg, endSeg do locss[i] = structure.GetSecondaryStructure( i ) end
    for i=startSeg, endSeg do
        if structure.GetSecondaryStructure( i ) == 'E' then 
            killedSS = true 
            structure.SetSecondaryStructure( i, 'L' )
        end 
    end

    selection.DeselectAll()
    save.Quicksave( QS_ShortUseTemp1 )   -- holds our "running" pose
    save.Quicksave( QS_ShortUseTemp2 )  -- holds our "best so far" result

    bestRescore = -10000.0
    for i=1, maxTries do
        local idxListInner = { }
        save.Quickload( QS_ShortUseTemp1 )  -- bring up our "running" pose

        selection.SelectRange( startSeg, endSeg )
        RebuildSelected( i ) -- not only up to maxTries times, but up to maxTries iterations per time
        selection.DeselectAll( )

        if math.abs( getScore( ) - startScore ) >= 0.01 then
              -- give the whole thing a little bit of shake-and-wiggle-time (deliberately not complete)
            structure.WiggleSelected( iters )
            getSegsWithinRangeOfTube( idxListInner, startSeg, endSeg, REBUILD_INNER_TUBE_RADIUS )
            ShakeOrMutate( iters/2, idxListInner )  -- only shake the "nearby" segments
            structure.WiggleAll( iters, true, true )  -- just a short bit of time on this wiggle to settle things down a bit

            rescore = getScore()
            if rescore > startScore then   -- we have a real improvement!
                save.Quicksave( QS_ShortUseTemp1 )  -- replace our "running" pose with this better one
                save.Quicksave( QS_ShortUseTemp2 ) -- our best-so-far is better
                startScore = rescore
                bestRescore = rescore
                changeSucceeded = true
            else    -- not improved, but maybe it's worth keeping anyway?
                if rescore > bestRescore and IsNeighborDifferentFromOld( oldScoreMat, 0.1 ) then
                    save.Quicksave( QS_ShortUseTemp2 )   -- our best-so-far is better
                    bestRescore = rescore
                    if startScore - rescore < maxFall then changeSucceeded = true end
                end
            end
        end
    end

    save.Quickload( QS_ShortUseTemp2 )
    if killedSS then
        for i=startSeg, endSeg do structure.SetSecondaryStructure( i, locss[ i ] ) end
    end
    selection.DeselectAll( )
    return changeSucceeded
end

function RebuildSection_Quaker( startSeg, endSeg, maxTries, maxFall )
    print( "RebuildSection_Quaker from "..startSeg.." to "..endSeg )
    local initScore = getScore( )
    local tries = 0
    repeat
        selection.DeselectAll( )
        selection.SelectRange( startSeg, endSeg )
        tries = tries + 1
        RebuildSelected( tries )
        newScore = getScore( )
    until newScore ~= initScore or tries > maxTries
    if newScore == initScore then 
        return false      -- we failed
    end

    recentbest.Save( )   -- N.B. This is exactly the sort of situation recentbest was made for!
    RebuildSelected( maxTries - 1 )   -- maybe use REBUILD_MAXITERS - 1 ?
    recentbest.Restore( )

    WiggleAndShakeForRebuilder( initScore, maxFall )
    -- quaking rebuild put a perform-quake-if-score-not-better-but-not-too-bad here
    return true
end

function RebuildSimple( idx, len, addCleanupWiggle )   -- randomizeR calls with idx=randomSeg(), len=random( 4 ), true
    save.Quicksave( QS_ShortUseTemp1 )
    local initscore = getScore( )

    selection.DeselectAll( )
    selection.SelectRange( idx, idx + len - 1 )
    RebuildSelected( 1 )

    structure.ShakeSidechainsSelected( 1 )
    structure.IdealizeSelected( )
    structure.ShakeSidechainsSelected( 1 )
    structure.WiggleAll( 25 )
    structure.ShakeSidechainsAll( 1 )
    structure.WiggleAll( 25 )
    local currscore = getScore( )
    if currscore > initscore then
        save.Quicksave( QS_ShortUseTemp1 )
        initscore = currscore
    end
    if addCleanupWiggle then 
        BlueFuse( )
        if currscore > initscore then
            save.Quicksave( QS_ShortUseTemp1 )
        end
    end
    save.Quickload( QS_ShortUseTemp1 )
    return true
end
----------------------------------------------------------------------
-- OPERATION-SPECIFIC CHOOSING FUNCTIONS
----------------------------------------------------------------------
function PickRebuildOperation( startSeg, endSeg, maxFall )
    local iters = QUICKWIGGLE_ITERS
    changeSucceeded = false

    if random( ) < 0.50 then
        algName = "RebuildInRange"
        local initScore = getScore( )
        changeSucceeded = RebuildInRange( startSeg, endSeg, iters, REBUILD_TRYCOUNT, maxFall )
        if changeSucceeded then
            if REBUILD_ADDIDEALIZE then
                DebugPrint("Putting Idealize after rebuild")
                save.Quicksave( QS_ShortUseTemp2 )  -- so we can throw out the idealize if we don't like result
                local score = getScore( )
                local check = IdealizeInRange( startSeg, endSeg, iters, REBUILD_USECUTS )
                if check and score < getScore( ) then 
                    save.Quicksave( QS_ShortUseTemp2 )
                    score = getScore( )
                end
                save.Quickload( QS_ShortUseTemp2 )
            end

            local idxListOuter = {}
            getSegsWithinRangeOfTube( idxListOuter, startSeg, endSeg, REBUILD_OUTER_TUBE_RADIUS )
            IterativeWiggleShakeWiggle( idxListOuter, iters  )
        end
    else
        -- this one comes with its own wiggler and skaker built-in
        algName = "RebuildSectionQ"      -- we'll just overrule alg's name choice for the time being
        changeSucceeded = RebuildSection_Quaker( startSeg, endSeg, 8, maxFall )  -- 8 => REBUILD_TRYCOUNT ?
    end

    return changeSucceeded
end

function PickPointRebuildOperation( idx, maxRange, maxFall )
    if maxFall == nil then maxFall = REBUILD_MAX_MEDIUM_SCOREDROP end
    if maxRange == nil then maxRange = math.floor( REBUILD_MAX_SEGLEN / 2 ) end
    local startIdx = math.max( 1, idx - random( maxRange ) )
    local endIdx = math.min( segCt, idx + random( maxRange ) )

    for i=idx, endIdx do
        if not IsMovableSeg(i) then 
            endIdx = i - 1
            break
        end
    end
    for i=idx, startIdx, -1 do
        if not IsMovableSeg(i) then 
            startIdx = i + 1
            break
        end
    end
    return PickRebuildOperation( startIdx, endIdx, maxFall )
end

function PickIdealizeRangeOperation( idx )
    local iters = QUICKWIGGLE_ITERS
    local startIdx, endIdx = GetRegionForOperation( idx )
    if endIdx - startIdx + 1 > IDEALIZE_MAX_SEGLEN then 
         if idx - startIdx + 1 > IDEALIZE_MAX_SEGLEN then 
             startIdx = math.max( 1, idx - IDEALIZE_MAX_SEGLEN) 
        end
    end
    if endIdx - startIdx + 1 > IDEALIZE_MAX_SEGLEN then 
         endIdx = startIdx + IDEALIZE_MAX_SEGLEN - 1
    end
    changeSucceeded = IdealizeInRange( startIdx, endIdx, iters, IDEALIZE_USE_CUTS )
    if changeSucceeded then qStabWiggle( iters, random() < PROB_QSTAB_EXTRA ) end

    return changeSucceeded
end

function PickIdealizePointOperation( idx )
    local iters = QUICKWIGGLE_ITERS
    changeSucceeded = IdealizeInRange( idx, idx, iters, IDEALIZE_USE_CUTS )
    if changeSucceeded then qStabWiggle( iters, random() < PROB_QSTAB_EXTRA ) end

    return changeSucceeded
end

--------------------------------------------------------------------------
-- MAIN AND SETUP AND CLEANING
--------------------------------------------------------------------------
function CleanPuzzleState( )
    SetCI( 1.0 )
    selection.DeselectAll()
    freeze.UnfreezeAll()
    DelBands()
    for i=1,segCt do
        if IsMovableSeg(i) then DeleteCutAtIndex( i ) end -- we can try, at least!
    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( " 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( )
    print("")
    ReportStats( )
    print( "TotalShakeTime="..TotalShakeTime )   -- because I want to know.
    CleanPuzzleState( )
end

function InitializePuzzleState()
    seedRandom( )
    InitialScore = getScore( )
    StartTime = os.time( )

    -- learn a few things
    FindAllMovableSegs()
    getHelices( )
    getSheets( )
    getLoops( )
    PuzzleAllowsMutate = IsPuzzleMutable( )
    if PuzzleAllowsMutate then IsDesignerPuzzle = checkAla() end
    CheckForContactMap()
    if PuzzleHasContactMap then InitializeContactMapSegList( CONTACTMAP_THRESHOLD ) end

    InitialClashImportance = behavior.GetClashImportance( )
    SCALE_MAXCI = InitialClashImportance

    selection.DeselectAll()
    freeze.UnfreezeAll()
    DelBands( )

    recentbest.Save( )
end

function main()
    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
        PerformOperation(  )
        didSomething = true
        PrintState( )
    end

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

xpcall( main, End )

Comments

[KarenCH]

some bug fixes from 2.0 - an attempt to regularize names before publishing had been incompletely done.

[KarenCH]

Notes for use:

1. It does better if the option of insisting on contact scores improving is turned off - the improvements that get rejected generally turn out to have found something worthwhile. I've taken to figuring that I'll get the contact points back later.

2. I usually have fusing turned off in early runs, and on in later ones.

3. Other than that, there doesn't seem to be much reason to change the default settings.