Recipe: Local Quake 2.0

created by spvincent

Profile

• Name: Local Quake 2.0
• ID: 34760
• Shared with: Public
• Parent: None
• Children:
  • Local Quake 2.0 504
  • Local Quake 2.0 for 504
  • Local Quake 2.0 start35
• Created on: November 04, 2011 at 00:23 AM UTC
• Updated on: November 04, 2011 at 00:23 AM UTC
• Description:

  Converted to V2 standard. Added some refinements.

Code

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

max_bands_per_residue = 6           -- Increase this for more bands at each segment
start_idx = 0                                         -- Set to a valid residue to start with that residue: otherwise 
                                                               -- start with a semi-random one,

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,

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


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

n_quakes  = 50000                         -- Go crazy. Effectively, quake forever.

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

monit = false

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


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

n_candidates = 0
best_score = 0

prev_band_strength = 0
prev_total_length = 0
prev_loss_after_banded_wiggle = 0
default_band_strength = 0.5

function roundto3 ( i )

    -- printing convenience

    if ( i == 0 ) then
        return 0
     else
        return i - i % 0.001
    end

end

function GetScore ()

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

end

function Coprime ( n )

    if ( n < 31 )  then
        return 1
    elseif  ( n >= 31*37 ) then
        return 1
    elseif ( n%31 ~= 0 ) then
        return 31
    else
        return 37
    end

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

    n_segments = 0
    segment_ids = {}
    segment_distances = {}

    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_length > 0 and total_band_length > 0 ) then
         if (  ( prev_total_length / total_band_length ) *  ( target_wiggle_reduction / prev_loss_after_banded_wiggle ) > 1 ) then
              band_strength = prev_band_strength * ( 1 + 0.15 * pseudorandom ( score_after_wiggle ) )
              
              if ( band_strength > 1 ) then
                   band_strength = 1
              end
           else
              band_strength = prev_band_strength * ( 1 - 0.15 *pseudorandom ( score_after_wiggle ) )
        end
      else
         band_strength = default_band_strength
    end

    if ( monit ) then
        print ( "bl " .. roundto3 ( total_band_length )  .. " pbl " .. roundto3 ( prev_total_length ) )
        print ( "bs " .. roundto3 ( band_strength ) .. " bps " .. roundto3 ( prev_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

   prev_band_strength = band_strength
   prev_total_length = total_band_length

end

function  Quake (  seg )

    selection.SelectAll ()

    init_score =  GetScore ()
    structure.WiggleSelected ( 1 )

    score_after_wiggle = GetScore ()

    prev_loss_after_banded_wiggle = init_score - score_after_wiggle

    if ( monit ) then
       print (  "seg ".. seg .. " score loss "  .. roundto3 ( prev_loss_after_banded_wiggle )  )
    end

    delete_all_bands ()

    structure.ShakeSidechainsSelected ( 1 )
    
    gas = GetScore () - score_after_wiggle
    if ( monit ) then
      print ( "      gas ".. gas )
    end

    structure.WiggleSelected ( 8 )

end

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

recentbest.Save ()

best_score = GetScore ()

multiplier = current.GetExplorationMultiplier ()
if ( multiplier < 1 ) then
  multiplier = 1
end

inc = Coprime ( n_residues )

print ( " init score ".. roundto3 ( 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 = 0

for i =  1, n_quakes  do

   if ( n_quakes_without_improvement >  time_for_a_change  ) then
          r = pseudorandom ( score )
	second_offset = 1 + r * ( n_residues - 2 )
         second_offset = second_offset - second_offset % 1
         n_quakes_without_improvement = 0
   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 " .. roundto3 ( score ) )
         n_quakes_without_improvement = 0
       else
          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 ( second_offset ~= 0 ) then
       print ( "i: " .. i .. " seg " .. idx .. "," .. second_idx )
     else
       print ( "i: " .. i .. " seg ".. idx )
    end

    create_bands ()

    Quake ( idx )

  end

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

end 

Comments

[marie_s]

Most of the time, don't destroy your protein