Recipe: tau-extremal optimization

created by CharaLilith

Profile

• Name: tau-extremal optimization
• ID: 105093
• Shared with: Public
• Parent: None
• Children:
  • tau-extremal optimization FORK
• Created on: November 15, 2021 at 11:46 AM UTC
• Updated on: November 15, 2021 at 11:46 AM UTC
• Description:

  In extremal optimization, the worst [segment] is changed randomly to see if it improves the fitness of a candidate solution. Tau-extremal optimization adds a bit of randomness to reduce the risk of getting stuck in local extrema. A tau of 0 makes every segment equally likely to be selected, and a higher tau will bias towards the worst segment more and more. A tau of around 1.37 gives an 80% chance of choosing from the worst 10 segments, and a bit over 25% chance of choosing the very worst.

Code

--[[
    a tau of 1.3 means there's about a 89% chance of selecting from the bottom 10 elements,
    and 27% of selecting the worst.  lower values make the distribution more uniform.
]]--
Behaviors = {"HidingImportance", "PackingImportance", "WigglePower", "BackboneHBondImportance", "SidechainHBondImportance",
                 "ClashImportance", "PairwiseImportance"}
for i = 1,7 do
    Behaviors[Behaviors[i]] = behavior["Get"..Behaviors[i]]()
    Behaviors[i] = nil
end

SaveSlot = nil

-- utility functions

function NOP()
    return
end

function NewSeed()
    math.randomseed(math.random(1, os.time()))
end

function RestoreSettings()
    for k,v in pairs(Behaviors) do
        behavior["Set"..k](v)
    end
    if SaveSlot == nil then
        print("Loaded recent best.")
        recentbest.Restore()
    else
        print("Loaded quicksave from slot "..SaveSlot..".")
        save.Quickload(SaveSlot)
    end
end

function Cleanup(errorMSG)
    if errorMSG == nil or string.find(errorMSG, "Cancelled") then
      print("Total score gained: "..EO.total_gain)
    else
      print(errorMSG)
    end
    RestoreSettings()
end

function SaveCurrent()
    --look for an empty quicksave spot, starting after the shortcut ones
    for i=9,100 do
        if save.QuicksaveEmpty(i) then
            save.Quicksave(i)
            return i
        end
    end
    --just overwrite slot 9, then
    save.Quicksave(9)
    return 9
end

function table.print(t)
    for k,v in pairs(t) do print(" " .. k .. ": "..v) end
end

function table.compare(t)
    return function(key1, key2)
        key1 = ""..key1
        key2 = ""..key2
        if t[key1] == nil then
          return key2
        end
        if t[key2] == nil then
          return key1
        end
        return t[key1] < t[key2]
    end
end


function neighbors(index)
    -- this actually returns the neighbors' neighbors, for freezing purposes
    local lower = nil
    local upper = nil
    local last = structure.GetCount()
    if index <= 2 then
        lower = last-(index-1)
    else
        lower = index -2
    end
    if index >= (last-1) then
        upper = 2-(last-index)
    else
        upper = index +2
    end
    return lower, upper
end



-- end utility functions



EO = {tau = 1.3, original = "current", worst_segments = {}, probabilities = {}, bands = {}, cuts = {}, freezes = {}, total_gain = 0.0, }

function EO:add_random(segment, options)
    -- use segment-1 and segment+1 as the other axes, and use random coordinates
    local theta = math.acos ( 2 * math.random () - 1 )
    local phi = 2 * math.pi * math.random ()
    local rho = options.max_band_length * math.random () ^ ( 1 / 3 ) + 0.001
    local lower,upper = neighbors(segment)
    if Options.cut_not_freeze == true then
        structure.InsertCut(lower)
        structure.InsertCut(upper)
        table.insert(self.cuts, lower)
        table.insert(self.cuts, upper)
    else
        freeze.Freeze(lower, true, false)
        freeze.Freeze(upper, true, false)
        table.insert(self.freezes, lower)
        table.insert(self.freezes, upper)
    end
    b = band.Add(segment, lower, upper, rho, theta, phi)
    if b == 0 then b = self:add_random(segment, options) end
    table.insert(self.bands, b)
    band.SetStrength(b, options.min_band_length + math.random() * (options.max_band_strength - options.min_band_length))
    band.SetGoalLength(b, options.min_band_goal_length + math.random()* (options.max_band_goal_length - options.min_band_goal_length))
    return b
end
function EO:remove_bands()
    count = band.GetCount()
    for k,b in pairs(self.bands) do
        if b < count then
            band.Delete(b)
        end
        table.remove(self.bands, k)
    end
end
function EO:remove_cuts()
    --assumes cuts can be healed
    for k,c in pairs(self.cuts) do
        structure.DeleteCut(c)
        table.remove(self.cuts, k)
    end
end
function EO:unfreeze()
    for k,f in pairs(self.freezes) do
        freeze.Unfreeze(f, true, false)
        table.remove(self.freezes, k)
    end
end

function EO:add_random_between_segments(s, options)
    NewSeed()
    local other = s
    while other == s do
        other = math.random(structure.GetCount())
    end
    b = band.AddBetweenSegments(s, other)
    if b == 0 then b = self:add_random_between_segments(s, options) end
    table.insert(self.bands, b)
    band.SetStrength(b, options.min_band_strength + math.random() * (options.max_band_strength - options.min_band_strength))
    band.SetGoalLength(b, options.min_band_goal_length + math.random() * (options.max_band_goal_length - options.min_band_strength))
    return b
end

function EO:calculate_probabilities()
    local sum = 0.0
    for i,v in ipairs(self.worst_segments) do
        if structure.IsLocked(i) then
            self.probabilities[""..v] = 0.0
        else
            self.probabilities[""..v] = (i+1.0) ^ -self.tau
        end
      sum = sum + self.probabilities[""..v]
    end
    return sum
end

function EO:chooseWithTau()
    --currently I'm recalculating these every time.  might (need to) improve that later.
    self.total_prob = self:calculate_probabilities()
    --table.print(probs)
    --print("total probability: "..total)
    NewSeed()
    local s = math.random()
    print("random threshold: "..s)
    for i=1,structure.GetCount(),1 do
      v = self.worst_segments[i]
      print("Testing segment "..v.." with P "..self.probabilities[""..v]/self.total_prob ..".")
      s = s - (self.probabilities[""..v] / self.total_prob)
      print("Remaining threshold: "..s)
      if s <= 0.0 then
        print("threshold met with "..v)
        return v
      end
    end
    return self.worst_segments[1]
end
    

function EO:run()
    local p = {} -- maps segment index to probability
    local p_score = current.GetScore()
    local b = nil
    self.worst_segments = {}

    for i=1,structure.GetCount() do
      p["" .. i] = current.GetSegmentEnergyScore(i)
      table.insert(self.worst_segments, i)
    end
    table.sort(self.worst_segments, table.compare(p))
    --worst = a[1]
    --print("Worst scoring segment: " ..worst.. ".  Score: " .. p[""..worst])
    
    local chosen = self:chooseWithTau()
    print("Segment "..chosen.." chosen with probability "..self.probabilities[""..chosen]/self.total_prob)
    b = -1
    if self.Options.band_segments then
        b = self:add_random_between_segments(chosen, self.Options)
    else
        b = self:add_random(chosen, self.Options)
    end
    if b ~= -1 then
        table.insert(self.bands, b)
    end
    if self.Options.set_wiggle_power then
        if behavior.HighPowerAllowed() then
            behavior.SetWigglePower("h")
        else
            behavior.SetWigglePower("m")
        end
    end
    behavior.SetClashImportance(0.05)
    behavior.SetSidechainHBondImportance(0.0)
    selection.Select(chosen)

    if self.Options.mutate and structure.IsMutable(chosen) then
        structure.MutateSidechainsSelected(chosen)
    end

    structure.WiggleAll(5)
    selection.Deselect(chosen)
    self:unfreeze()
    self:remove_cuts()
    --for i=1,structure.GetCount() do
    --    structure.SetSecondaryStructure(i, "a")
    --end
    self:remove_bands()
    structure.ShakeSidechainsAll(1)
    
    if self.Options.mutate_all then
        structure.MutateSidechainsAll(1)
    end

    behavior.SetPackingImportance(1.4)
    behavior.SetHidingImportance(1.8)
    behavior.SetClashImportance(1.0)
    behavior.SetSidechainHBondImportance(0.7)
    structure.WiggleAll(5)
    structure.ShakeSidechainsAll(2)
    behavior.SetPackingImportance(1.0)
    behavior.SetHidingImportance(0.8)
    behavior.SetSidechainHBondImportance(1.5)
    if self.Options.set_wiggle_power then
        behavior.SetWigglePower("a")
    end
    structure.WiggleAll(5)
    print("Segment score changed from "..p[""..chosen].." to " ..current.GetSegmentEnergyScore(chosen))

    local gain = current.GetScore() - p_score
    print("Score gain of "..gain)
    if gain < 0.0 then
        save.Quickload(SaveSlot)
    else
        save.Quicksave(SaveSlot)
        print("Positive score gain!  New solution saved to quicksave slot "..SaveSlot)
        self.total_gain = math.max(self.total_gain, gain)
    end
    --remove only added bands
    self:remove_bands()
    return math.max(p_score, gain+p_score)
end
    
function main()
    Options = dialog.CreateDialog("tau-extremal optimization options")
    Options.tau_explanation = dialog.AddLabel("Values between 1.1 and 1.7 are recommended.")
    Options.tau = dialog.AddSlider("tau", 1.3, 0.0, 2.0, 2)
    Options.R_label = dialog.AddLabel("Length of space bands")
    Options.min_band_length = dialog.AddSlider("min", 1.0, 0.0, 100.0, 1)
    Options.max_band_length = dialog.AddSlider("max", 10.0, 1.0, 100.0, 1)
    Options.goal_label = dialog.AddLabel("Goal length of rubber bands")
    Options.min_band_goal_length = dialog.AddSlider("min", 0.0, 0.0, 100.0, 1)
    Options.max_band_goal_length = dialog.AddSlider("max", 10.0,0.0, 100.0, 1)
    Options.strength_label = dialog.AddLabel("Strength of rubber bands")
    Options.min_band_strength = dialog.AddSlider("min", 0.0, 0.0, 10.0, 1)
    Options.max_band_strength = dialog.AddSlider("max", 5.0, 0.0, 10.0, 1)
    Options.band_segments = dialog.AddCheckbox("Band segments together?", false)
    Options.mutate = dialog.AddCheckbox("Mutate worst segment?", false)
    Options.mutate_all = dialog.AddCheckbox("Mutate all segments?", false)
    Options.iterations = dialog.AddSlider("Iterations", 1000, 1, 10000, 0)
    Options.explain_score = dialog.AddLabel("If you check this, then iterations are ignored.")
    Options.stop_after_score = dialog.AddCheckbox("Stop after minimum score gain", false)
    Options.minimum_score_gain = dialog.AddTextbox("Minimum score gain", 1000.0)
    Options.set_wiggle_power = dialog.AddCheckbox("Adjust wiggle power?", true)
    Options.cut_not_freeze = dialog.AddCheckbox("Cut instead of freeze?", false)
    Options.OK = dialog.AddButton("OK", 1)
    Options.Cancel = dialog.AddButton("Cancel", 0)
    if dialog.Show(Options) <= 0 then
        return Cleanup("User canceled.")
    end
    EO.Options = {}
    for i,v in ipairs(Options._Order) do
        if Options[v].controlType ~= 0 then
            if Options[v].value ~= nil then
                print(v .. ": "..tostring(Options[v].value))
            end
            EO.Options[v] = Options[v].value
        end
    end

    behavior.SetHidingImportance(2.0)
    behavior.SetPackingImportance(0.7)

    SaveSlot = SaveCurrent()

    if EO.Options.stop_after_score then
        while EO.total_gain < tonumber(EO.Options.minimum_score_gain) do
            EO:run()
        end
    else
        for i=1,EO.Options.iterations do
            EO:run()
        end
    end
    
    Cleanup()
end
xpcall(main, Cleanup)

Comments

[zo3xiaJonWeinberg]

Nothing to do with half of pi, right? But works

[CharaLilith]

Before selecting "Stop after minimum score gain" stopped counterintuitively. This is because I was tracking the total score gained rather than just the maximum. It should be easier to set it and forget it now.

[CharaLilith]

okay, now there's a bug where it crashes on line 119 (band.add_random_between_segments) saying that local 'options' is nil. I cannot figure out why options would be nil, or why it wouldn't happen for a while. Maybe it gets GCed for some reason?

anyway, that won't get fixed unless 1. I suddenly become way smarter, or 2. someone else figures out the issue.

[alcor29]

Can this be modified to not choose locked segments, as for example in binder problems?

[LociOiling]

Sorry for the slow reply, but here goes.

The function band.add_random_between_segments is defined starting at line 111.

It has two arguments.

Then it calls itself at line 118, passing only one argument.

So I'd suggest try changing line 118:

if b == 0 then b = band.add_random_between_segments(s) end

to:

 
if b == 0 then b = band.add_random_between_segments(s,options) end

Haven't tried it yet, just a thought so far.

[LociOiling]

I tested the suggested fix on puzzle 2070. I had "band segments together" checked, along with both mutate options. The recipe ran this way for over two hours with no crashes, and it found some points.

[CharaLilith]

I will look into that, but probably.

[CharaLilith]

Thanks for the help!