Recipe: Local Quake 3.1

created by spvincent

Profile

• Name: Local Quake 3.1
• ID: 44398
• Shared with: Public
• Parent: None
• Children: None
• Created on: October 21, 2012 at 23:48 PM UTC
• Updated on: October 21, 2012 at 23:48 PM UTC
• Description:

  See first comment

Code

-- Bands radiate from two segments to spatially close neighbours
-- Foregoes Fuses and generally tries to optimize for speed

-- Some default values, modifiable by the initial dialog

max_bands_per_residue = 6           -- Increase this for more bands at each segment

second_offset = 3                              -- The difference between the banding centres.
                                                               -- Set to 0 for a single centre. This is the most likely one to change but it'll
                                                               -- change anyway after time_for_a_change (see below) quakes without improvement.

target_wiggle_reduction = 150        -- On wiggle, try and reduce the score by this amount after creating the bands. 
                                                                -- If it's wildly out, reduce/increase band strength
                                                                -- for the next try (don't try to adjust for the current scramble).

max_segment_length = 20               -- Don't create bands greater than this length. Set to 999 or something if you're 
                                                               -- not worried about the potential for long bands creeping in,

mutate_not_shake = false                -- for design puzzles

-- Other parameters, less likely to be changed and so not added in the dialog to avoid clutter

time_for_a_change = 10                  -- After doing this number of quakes without any improvement, change the value of 
                                                               -- second_offset to something random.

scale_band_strength = true            -- If true, band strength will be scaled inversely with length

max_allowable_band_strength = 3

---          Some parameters you'll probably not want to change

use_ranked_score = 1                       -- if 1, use the score being used for ranking (best normally). If 0, use the
                                                               --  old style scoring system

band_to_minima = true                    -- Band to minima or maxima in the distance. See bands_from_segment for meaning. 

min_exclusion =  999                       -- These can be set to establish an "exclusion zone", within which no bands will
                                                             --  be created. This might be used to prevent banding to a wayward loop.
                                                              --  If min_exclusion > max_exclusion no exclusion will be set.
max_exclusion =  1

start_idx = 0                                        -- Set to a valid residue to start with that residue: otherwise 
                                                              -- start with a semi-random one,

monit = false                                        -- prints out performance info

---------- A few more globals


candidate_ids_start = {}
candidate_ids_end = {}
candidate_distances = {}

n_candidates = 0
best_score = 0

gain_after_shake = 0

band_strength = 0
prev_band_strength = 0

total_band_length = 0
prev_total_band_length = 0

loss_after_banded_wiggle = 0
prev_loss_after_banded_wiggle = 0

default_band_strength = 0.3
wiggle_accuracy = 3.0

function r3 ( x )

-- Round to 3 decimal places

  t  = 10 ^ 3

  return math.floor ( x*t + 0.5 ) /  t 

end


function GetScore ()

  if ( use_ranked_score == 1 ) then
     return current.GetScore ()
   else
     return current.GetEnergyScore ()
  end

end

function Coprime ( 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}
    
    i = #primes
   -- find the highest prime < 40% of the residue count and which is coprime with the number of residues

    while ( i >= 1 ) do
        if ( primes [ i ] < n*0.4 and primes [ i ] % n ~= 0 ) then
           return primes [ i ] 
       end

       i = i - 1

    end

    return 1

end

function pseudorandom ( score )

    -- Returns a pseudorandom number >= 0 and < 1.
    -- Based on the fractional part of the score

    if ( score >= 0 ) then
      return  score % 1
     else
      return (-score ) % 1
    end

end

function  delete_all_bands ()

    band.DeleteAll ()

 end

function IsPuzzleMutable ()

    for i = 1 , n_residues do
       if ( structure.IsMutable ( i ) ) then
           return true
       end
    end

    return false

end

function  ShellSort ( ids , distances , n )

    -- Adapted from Numerical Recipes in C

    local  inc = 1

    repeat
      inc = inc * 3 + 1
    until inc > n

    repeat

        inc = inc / 3
        inc = inc - inc % 1

        for i = inc + 1 , n do

            v = distances [ i ]
            w = ids [  i ]

            j = i
            flag = false

            while ( flag == false and distances [ j - inc ]  > v ) do

                distances [ j ] = distances [ j - inc ]
                ids [ j ] = ids [ j - inc ]

                j = j - inc

                if ( j <= inc ) then
                  flag = true
                end

            end

          distances [ j ] = v
           ids [ j ] = w

        end

    until inc <= 1

end

function test_for_minimum ( i , seg , seg_distances )

   if (  ( i ~= seg ) and
         ( seg_distances [ i ] <= seg_distances [ i - 1 ]  ) and 
         ( seg_distances [ i ] <= seg_distances  [ i - 2 ] ) and
         ( seg_distances [ i ] <= seg_distances [ i + 1 ] ) and
         ( seg_distances [ i ] <= seg_distances [ i + 2 ] ) ) then
     return true
    else
     return false
  end
end

function test_for_maximum ( i , seg , seg_distances )

   if ( ( seg_distances [ i ]  >= seg_distances [ i - 1 ]  ) and 
         ( seg_distances [ i ] >= seg_distances  [ i - 2 ] ) and
         ( seg_distances [ i ] >= seg_distances [ i + 1 ] ) and
         ( seg_distances [ i ] >= seg_distances [ i + 2 ] ) ) then
     return true
    else
     return false
  end

end

function is_segment_moveable ( seg )

    seg_backbone_frozen, seg_sidechain_frozen = freeze.IsFrozen ( seg )

    if ( seg_backbone_frozen == true ) then
      return false
    end

    s = structure.IsLocked ( seg )
    
    if ( s == true ) then 
      return false     -- greyed out as in a design puzzle
    else
      return true
    end

end

function get_band_segment_data ( seg )

    seg_distances = {}

    -- Get the function that looks like this: 
    -- x = residue number
    -- y = distance of that segment from seg
    -- ignore the degenerate minimum

    -- Then look for minima in that function as candidates for banding.

    n_segments = 0
    segment_ids = {}
    segment_distances = {}

    seg_ok = is_segment_moveable ( seg )
    if ( seg_ok == false ) then
       return
    end

    for i = 1, n_residues do
        segment_distances [ i ] =  structure.GetDistance  ( seg , i )
    end

    for  i = 3 , n_residues - 3  do
        
       k = segment_distances [ i ]

        if ( ( band_to_minima == true ) and 
              ( test_for_minimum ( i , seg , segment_distances ) == true ) ) then
                  n_segments = n_segments + 1
                  segment_ids  [ n_segments ]  = i
                  segment_distances [ n_segments ] =  segment_distances [ i ]
          elseif ( ( band_to_minima == false ) and 
                        ( test_for_maximum ( i , seg , segment_distances ) == true )  ) then
                            n_segments = n_segments + 1
                           segment_ids  [ n_segments ]  =  i
                           segment_distances [ n_segments ] =  segment_distances [ i ]
        end

    end

    ShellSort ( segment_ids , segment_distances , n_segments )

    if ( n_segments >= max_bands_per_residue ) then
      n_c = max_bands_per_residue
    else
      n_c = n_segments
    end

    for i = 1 , n_c do

        if ( ( segment_distances [ i ] <= max_segment_length ) and 
              ( seg < min_exclusion or seg > max_exclusion ) and
              ( segment_ids [ i ] < min_exclusion or segment_ids [ i ] > max_exclusion ) ) then
          n_candidates = n_candidates + 1
          candidate_ids_start [ n_candidates ] = seg
          candidate_ids_end [ n_candidates ] = segment_ids [ i ]
          candidate_distances [ n_candidates ] = segment_distances [ i ]
       end

    end

end

function create_bands ()

    total_band_length = 0
    for i = 1 , n_candidates do
         total_band_length = total_band_length + candidate_distances [ i ]
    end

    if ( scale_band_strength == true ) then
        average_band_length = total_band_length / n_candidates
    end

    -- Adaptively change band strength based principally on how much the previous wiggle overshot/undershot the target.

    if ( prev_total_band_length > 0 and total_band_length > 0 ) then
         if (  ( target_wiggle_reduction / prev_loss_after_banded_wiggle ) > 1 ) then
              band_strength = prev_band_strength * ( 1 + 0.15 * pseudorandom ( score_after_banded_wiggle ) )
              
              band_strength = math.min ( band_strength , max_allowable_band_strength )

           else
              band_strength = prev_band_strength * ( 1 - 0.15 *pseudorandom ( score_after_banded_wiggle ) )
        end
      else
         band_strength = default_band_strength
    end

    for i = 1 , n_candidates do

        band.AddBetweenSegments ( candidate_ids_start [ i ] , candidate_ids_end [ i ] )

        if ( scale_band_strength == true ) then
             band.SetStrength ( i , band_strength * average_band_length /  candidate_distances [ i ] )
          else
             band.SetStrength ( i , band_strength )
        end

    end

end

function  Quake (  seg )

   selection.SelectAll ()

    init_score =  GetScore ()

    -- behavior.SetWiggleAccuracy ( 2.0 )
    structure.WiggleSelected ( 1 )

    score_after_banded_wiggle = GetScore ()

    loss_after_banded_wiggle = init_score - score_after_banded_wiggle

    delete_all_bands ()

    if ( loss_after_banded_wiggle < 8* target_wiggle_reduction ) then
    
        if ( mutate_not_shake == true ) then
            structure.MutateSidechainsSelected ( 1 )
          else
            structure.ShakeSidechainsSelected ( 1 )
        end
   
        gain_after_shake = GetScore () - score_after_banded_wiggle

        -- behavior.SetWiggleAccuracy ( wiggle_accuracy )
        structure.WiggleSelected ( 8 )

    end

end

function get_parameters ()

  local dlog = dialog.CreateDialog ( "Local Quake 3.1" )

  dlog.target_reduction = dialog.AddSlider ( "Target reduction" , target_wiggle_reduction , 10 , 1000 , 1 )  
  dlog.n_bands = dialog.AddSlider ( "Max bands" , max_bands_per_residue , 1 , 12 , 0 )  
  dlog.seg_length = dialog.AddSlider ( "Max seg length" , max_segment_length , 1 , 60 , 1 )  

  -- dialog.AddLabel ( "--------------------" )

  dlog.seg_offset = dialog.AddSlider ( "Second offset" , second_offset , 0 , n_residues - 1  , 0 )  
  dlog.change_after = dialog.AddSlider ( "Change after" ,  time_for_a_change , 1 , 40  , 0 )  
  dlog.inv_scale= dialog.AddCheckbox ( "Inverse scale" , scale_band_strength )  

  design_puzzle = IsPuzzleMutable ()

  if ( design_puzzle == true ) then
     dlog.mutate= dialog.AddCheckbox ( "Mutate not shake" , false )  
  end
    
  dlog.ok = dialog.AddButton ( "OK" , 1 )
  dlog.cancel = dialog.AddButton ( "Cancel" , 0 )

  if ( dialog.Show ( dlog ) > 0 ) then
    
    target_wiggle_reduction = dlog.target_reduction.value
    max_bands_per_residue = dlog.n_bands.value 
    max_segment_length = dlog.seg_length.value 

    second_offset = dlog.seg_offset.value
    time_for_a_change = dlog.change_after.value
    scale_band_strength = dlog.inv_scale.value
    
    if ( design_puzzle == true ) then
       mutate_not_shake = dlog.mutate.value
    end

    return true
  else
    return false
 end

end

function main ()

    n_residues = structure.GetCount ()
    band_strength = 0.5
    behavior.SetClashImportance ( 1 )

    -- wiggle_accuracy = behavior.GetWiggleAccuracy ()

    if ( get_parameters () == false ) then
      return               -- graceful exit
    end

    recentbest.Save ()

    best_score = GetScore ()

    multiplier = math.max ( current.GetExplorationMultiplier () , 1 )

    inc = Coprime ( n_residues )

    print ( " init score ".. r3 ( best_score ) )

    if ( start_idx < 1 or start_idx  > n_residues ) then
        r = pseudorandom ( best_score )
        start_idx = 1 + n_residues * r
        start_idx = start_idx - start_idx % 1
    end

    idx = start_idx
    n_quakes_without_improvement = 1
    score_after_quake = 0
    n_candidates = 0

    i = 1

    while true do

        if ( n_quakes_without_improvement >=  time_for_a_change  ) then
             r = pseudorandom ( score_after_quake )
             second_offset = 1 + r * ( n_residues - 2 )
             second_offset = second_offset - second_offset % 1
             n_quakes_without_improvement = 1
        end

        recentbest.Restore  ()
  
         -- Need to record any improvements from the previous iteration

        score = GetScore ()
        if ( score > best_score ) then
            best_score = score
            print ( "Improvement to " .. r3 ( score ) )
            n_quakes_without_improvement = 0
          elseif ( n_candidates > 0 ) then
            n_quakes_without_improvement = n_quakes_without_improvement + 1
        end

        delete_all_bands ()

        n_candidates = 0
        get_band_segment_data  ( idx )
    
        if ( second_offset ~= 0 ) then
            second_idx = idx + second_offset
            while ( second_idx > n_residues ) do
               second_idx = second_idx - n_residues
            end
    
            get_band_segment_data  ( second_idx )
  
        end

        if ( n_candidates > 0 ) then
    
            if ( monit == false ) then
                if ( second_offset ~= 0 ) then
                     print ( "i: " .. i .. " seg " .. idx .. "," .. second_idx )
                  else
                     print ( "i: " .. i .. " seg ".. idx )
                end
            end

            create_bands ()

            Quake ( idx )

            score_after_quake = GetScore ()

            if ( monit ) then
                print (  idx .. " , " ..
                            second_offset .. " , " .. 
                            r3 ( loss_after_banded_wiggle ) .. " , " ..
                            r3 ( gain_after_shake ) .. " , " ..
                            r3 ( band_strength ) .. " , " ..
                            r3 ( total_band_length ) )
            end

           prev_band_strength = band_strength
           prev_total_band_length = total_band_length
           prev_loss_after_banded_wiggle = loss_after_banded_wiggle

        end

        idx = idx  + inc
        if ( idx > n_residues ) then
            idx = idx - n_residues
        end
    
        i = i + 1

    end

end

function cleanup ()
   
    print ( "Cleaning up" )
    behavior.SetClashImportance ( 1.0 )
    recentbest.Restore  ()
 
end

xpcall ( main , cleanup )

Comments

[spvincent]

A couple of changes to make it behave a little better on design puzzles (Thanks Marie for pointing out the problem).