Code
-- Bands radiate from nodes 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
n_nodes = 2
distance_between_nodes = 3
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).
kMaxIterations = 10 -- Sanity check for the new wiggle loop
min_residue = 1 -- Change these if you want to rebuild part of a protein only. For example, to
-- rebuild 53-64, set min_residue to 53 and max_residue to 64. To restore, set to
-- 1 and 9999 (or any large number) respectively
max_residue = 999
max_band_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,
min_band_length = 2 -- Don't create bands less than this length.
mutate = false -- for design puzzles
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
kInitBandStrength = 1.0
kMaxAllowableBandStrength = 3
kSlackTolerance = 1.8 -- The idea is to compress the protein until a certain target reduction is met. In practice it's -- expedient to allow quite a bit of "wiggle room" in this value and kSlackTolerance defines -- how much. See code for usage but high values allow for more wiggle room.
--- Some parameters you'll probably not want to change
band_to_minima = true -- Band to minima or maxima in the distance. See bands_from_segment for meaning.
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
test_for_credit = false
---------- A few more globals
candidate_ids_start = {}
candidate_ids_end = {}
candidate_distances = {}
n_candidate_bands = 0
best_score = 0
gain_from_wiggle_or_mutate = 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
score_after_banded_wiggle = 0
default_band_strength = 0.3
---------- Helix data
helix_starts = {}
helix_ends = {}
use_helix = {}
n_helices = 0
kHelixThreshold = 4 -- doesn't matter if you mangle helices of this length or less
kMaxAllowableAngle = 30 -- max allowable angle in degrees from perpendicular
disallow_bands_parallel_to_helix = true
kTemp = 5 -- quicksave slots
function r3 ( x )
-- Round to 3 decimal places
t = 10 ^ 3
return math.floor ( x*t + 0.5 ) / t
end
function GetScore ()
return current.GetEnergyScore ()
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 n % primes [ i ] ~= 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 IsPuzzleMutable ()
for i = 1 , n_residues do
if ( structure.IsMutable ( i ) ) then
return true
end
end
return false
end
function CanAcceptImprovement ()
-- Done this way to work round a crashing bug in
-- creditbest.AreConditionsMet ()
if ( test_for_credit == false ) then
return true
elseif ( creditbest.AreConditionsMet () == true ) then
return true
else
return false
end
end
function GetHelices ()
within_helix = false
for i = 1, n_residues do
if ( structure.GetSecondaryStructure ( i ) == "H" ) then
if ( within_helix == false ) then
-- start of a new helix
within_helix = true
n_helices = n_helices + 1
helix_starts [ n_helices ] = i
end
elseif ( within_helix == true ) then
-- end of a helix
within_helix = false
helix_ends [ n_helices ] = i -1
end
end -- for i
if ( within_helix == true ) then
helix_ends [ n_helices ] = n_residues
end
end
function TestBandAgainstHelices ( band_a , band_b )
-- Tests the putative band to see, should band_a lie on a helix, whether the band is reasonably perpendicular to the helix.
if ( disallow_bands_parallel_to_helix == false ) then
return true
end
end_in_helix = false
local i = 1
while ( ( i <= n_helices ) and ( end_in_helix == false ) ) do
if ( ( band_a >= helix_starts [ i ] ) and
( band_a <= helix_ends [ i ] ) and
( use_helix [ i ] == true ) ) then
end_in_helix = true
-- Compute an approximate angle with the helix based on the cosine rule
if ( ( band_a - helix_starts [ i ] ) > ( helix_ends [ i ] - band_a ) ) then
-- The point lies closest to the helix end. Compare against the start of the helix.
a = structure.GetDistance ( band_a , helix_starts [ i ] )
b = structure.GetDistance ( band_a , band_b )
c = structure.GetDistance ( helix_starts [ i ] , band_b )
else
-- The point lies closest to the helix start. Compare against the end of the helix.
a = structure.GetDistance ( band_a , helix_ends [ i ] )
b = structure.GetDistance ( band_a , band_b )
c = structure.GetDistance ( helix_ends [ i ] , band_b )
end
cos_angle = ( a*a + b*b - c*c ) / ( 2*a*b )
-- We're looking for an angle of about 90 degrees give or take a few: other angles will be rejected
angle = math.deg ( math.acos ( cos_angle ) )
if ( ( angle < 90 - kMaxAllowableAngle ) or
( angle > 90 + kMaxAllowableAngle ) ) then
return false
else
return true
end
end
i = i + 1
end
return true
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 GetBandSegmentData ( 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_band_length ) and
( segment_distances [ i ] >= min_band_length ) and
( TestBandAgainstHelices ( seg , segment_ids [ i ] ) == true ) and
( TestBandAgainstHelices ( segment_ids [ i ] , seg ) == true ) ) then
n_candidate_bands = n_candidate_bands + 1
candidate_ids_start [ n_candidate_bands ] = seg
candidate_ids_end [ n_candidate_bands ] = segment_ids [ i ]
candidate_distances [ n_candidate_bands ] = segment_distances [ i ]
end
end
end
function CreateBands ()
total_band_length = 0
for i = 1 , n_candidate_bands do
total_band_length = total_band_length + candidate_distances [ i ]
end
if ( scale_band_strength == true ) then
average_band_length = total_band_length / n_candidate_bands
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.2 * pseudorandom ( score_after_banded_wiggle ) )
band_strength = math.min ( band_strength , kMaxAllowableBandStrength )
else
band_strength = prev_band_strength * ( 1 - 0.2 *pseudorandom ( score_after_banded_wiggle ) )
end
else
band_strength = default_band_strength
end
for i = 1 , n_candidate_bands 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
if ( i == 1 ) then
undo.SetUndo ( true )
end
end
undo.SetUndo ( true )
end
function Wiggle ( step , threshold )
n_it = 0
new_score = GetScore ()
repeat
prev_score = new_score
structure.WiggleSelected ( step )
new_score = GetScore ()
n_it = n_it + 1
until ( ( math.abs ( new_score - prev_score ) < threshold ) or ( n_it > kMaxIterations ) )
end
function AdjustBandStrength ( factor )
n_bands = band.GetCount ()
for i = 1 , n_bands do
strength = band.GetStrength ( i )
band.SetStrength ( i , strength * factor )
strength = band.GetStrength ( i )
if ( i == 1 ) then
undo.SetUndo ( false )
end
end
undo.SetUndo ( true )
end
function WiggleBandedProtein ()
-- We are given a protein with a bunch of bands between segments: the intent is to wiggle until
-- until the reduction from the start position is reasonably close to the desired value.
-- Wiggle once to get an approximate idea of whether bands need to be strengthened, weakened,
-- or whether they're just right.
step = 2
init_score = GetScore ()
structure.WiggleSelected ( step )
score = GetScore ()
loss_after_banded_wiggle = init_score - score
--print ( "WBP : init reduction " .. r3 ( loss_after_banded_wiggle ) )
if ( loss_after_banded_wiggle < 0 ) then
return 0 -- Rare
elseif ( loss_after_banded_wiggle > kSlackTolerance * target_wiggle_reduction ) then
-- Reduce band strength
return 1
elseif ( loss_after_banded_wiggle < target_wiggle_reduction / kSlackTolerance ) then
-- Increase band strength
return -1
else
return 0
end
end
function Compress ()
-- After creating the band , attempt to do the wiggle series (twice if necessary)
save.Quicksave ( kTemp )
s = WiggleBandedProtein ()
if ( s == 1 ) then
-- Bands too strong
save.Quickload ( kTemp )
multiplier = ( 1 - 0.4 * math.random () )
band_strength = band_strength * multiplier
AdjustBandStrength ( multiplier )
s = WiggleBandedProtein ()
elseif ( s == -1 ) then
-- Bands too weak
save.Quickload ( kTemp )
multiplier = ( 1 + 0.4 * math.random () )
band_strength = band_strength * multiplier
AdjustBandStrength ( multiplier )
s = WiggleBandedProtein ()
end
end
function Quake ( seg )
selection.SelectAll ()
Compress ()
band.DeleteAll ()
if ( mutate == true ) then
structure.MutateSidechainsSelected ( 1 )
else
structure.WiggleSelected ( 1 , false , true )
end
gain_from_wiggle_or_mutate = GetScore () - score_after_banded_wiggle
if ( ( mutate == true ) and ( math.abs ( gain_from_wiggle_or_mutate ) < 1e-5 ) ) then
-- Mutate didn't work: try wiggling sidechains
structure.WiggleSelected ( 1 , false , true )
end
Wiggle ( 4 , 0.1 )
end
function GetParameters ()
local dlog = dialog.CreateDialog ( "Local Quake 5.0" )
dlog.min_residue = dialog.AddSlider ( "Min residue" , 1 , 1 , n_residues , 0 )
dlog.max_residue = dialog.AddSlider ( "Max residue" , n_residues , 1 , n_residues , 0 )
dlog.target_reduction = dialog.AddSlider ( "Target reduction" , target_wiggle_reduction , 1 , 1000 , 0 )
dlog.n_bands = dialog.AddSlider ( "Max bands" , max_bands_per_residue , 1 , 12 , 0 )
dlog.min_seg_length = dialog.AddSlider ( "Min band length" , min_band_length , 1 , 80 , 1 )
dlog.max_seg_length = dialog.AddSlider ( "Max band length" , max_band_length , 1 , 80 , 1 )
-- dialog.AddLabel ( "--------------------" )
dlog.n_nodes = dialog.AddSlider ( "N nodes" , n_nodes , 1 , 4 , 0 )
dlog.distance_between_nodes = dialog.AddSlider ( "Distance between nodes" , distance_between_nodes , 1 , n_residues/2 , 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 )
dlog.helix_bands = dialog.AddCheckbox ( "Disallow helix ll bands" , disallow_bands_parallel_to_helix )
design_puzzle = IsPuzzleMutable ()
if ( design_puzzle == true ) then
dlog.mutate= dialog.AddCheckbox ( "Mutate" , false )
end
dlog.credit_check = dialog.AddCheckbox ( "Check conditions met " , test_for_credit )
dlog.ok = dialog.AddButton ( "OK" , 1 )
dlog.cancel = dialog.AddButton ( "Cancel" , 0 )
if ( dialog.Show ( dlog ) > 0 ) then
min_residue = dlog.min_residue.value
max_residue = dlog.max_residue.value
target_wiggle_reduction = dlog.target_reduction.value
max_bands_per_residue = dlog.n_bands.value
min_band_length = dlog.min_seg_length.value
max_band_length = dlog.max_seg_length.value
n_nodes = dlog.n_nodes.value
distance_between_nodes = dlog.distance_between_nodes.value
time_for_a_change = dlog.change_after.value
scale_band_strength = dlog.inv_scale.value
disallow_bands_parallel_to_helix = dlog.helix_bands.value
if ( design_puzzle == true ) then
mutate = dlog.mutate.value
end
test_for_credit = dlog.credit_check.value
return true
else
return false
end
end
function main ()
n_residues = structure.GetCount ()
band_strength = kInitBandStrength
behavior.SetClashImportance ( 1 )
if ( GetParameters () == false ) then
return -- graceful exit
end
print ( "Local Quake 5.0 : New Chapters" )
print ( "Nodes from " .. min_residue .. " to " .. max_residue )
print ( "Target reduction " .. target_wiggle_reduction )
print ( "Maximum of " .. max_bands_per_residue .. " bands with length between " .. min_band_length .. " and " .. max_band_length )
print ( n_nodes .. " nodes initially spaced " .. distance_between_nodes .. " segments apart" )
if ( mutate == true ) then
print ( "Mutating" )
end
recentbest.Save ()
best_score = GetScore ()
GetHelices ()
for i = 1 , n_helices do
if ( helix_ends [ i ] - helix_starts [ i ] > kHelixThreshold ) then
use_helix [ i ] = true
else
use_helix [ i ] = false
end
end
if ( max_band_length < min_band_length ) then
print ( "Swapping min/max band lengths" )
temp = max_band_length
max_band_length = min_band_length
min_band_length = temp
end
if ( max_residue < min_residue ) then
print ( "Rebuild range error." )
print ( "Max residue (" .. max_residue .. ") must be greater than or equal to Min residue (" .. min_residue .. ")" )
return
end
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_candidate_bands = 0
local nq = 1 -- printing purposes only
i = 1
while true do
if ( ( n_quakes_without_improvement >= time_for_a_change ) and ( n_nodes > 1 ) ) then
r = pseudorandom ( score_after_quake )
distance_between_nodes = 1 + r * ( n_residues - 2 )
distance_between_nodes = distance_between_nodes - distance_between_nodes % 1
n_quakes_without_improvement = 1
print ( "Node distance adjusted to " .. distance_between_nodes )
end
recentbest.Restore ()
-- Need to record any improvements from the previous iteration
score = GetScore ()
if ( ( score > best_score ) and
( test_for_credit == false or creditbest.AreConditionsMet () == true ) ) then
best_score = score
print ( "Improvement to " .. r3 ( score ) )
n_quakes_without_improvement = 0
elseif ( n_candidate_bands > 0 ) then
n_quakes_without_improvement = n_quakes_without_improvement + 1
end
band.DeleteAll ()
node_list = {}
n_candidate_bands = 0
for j = 1 , n_nodes do
node = idx + ( j -1 ) *distance_between_nodes
while ( node > n_residues ) do
node = node - n_residues
end
if ( node >= min_residue and node <= max_residue ) then
GetBandSegmentData ( node )
table.insert ( node_list , node )
end
end
if ( n_candidate_bands > 0 ) then
-- There has to be a better way of doing this
if ( #node_list == 1 ) then
print ( "i: " .. nq .. " (".. node_list [ 1 ] .. ")" )
elseif ( #node_list == 2 ) then
print ( "i: " .. nq .. " (".. node_list [ 1 ] .. "+" .. node_list [ 2 ] .. ")" )
elseif ( #node_list == 3 ) then
print ( "i: " .. nq .. " (".. node_list [ 1 ] .. "+" .. node_list [ 2 ] .. "+" .. node_list [ 3 ] .. ")" )
elseif ( #node_list == 4 ) then
print ( "i: " .. nq .. " (".. node_list [ 1 ] .. "+" .. node_list [ 2 ] .. "+" .. node_list [ 3 ] .. "+" .. node_list [ 4 ] ..")" )
elseif ( #node_list > 4 ) then
print ( "i: " .. nq .. " (".. node_list [ 1 ] .. "+" .. node_list [ 2 ] .. "+" .. node_list [ 3 ] .. "+" .. node_list [ 4 ] .. "+" ..")" )
end
nq = nq + 1
CreateBands ()
Quake ( idx )
score_after_quake = GetScore ()
if ( monit ) then
print ( " " ..
r3 ( loss_after_banded_wiggle ) .. " , " ..
r3 ( gain_from_wiggle_or_mutate ) .. " , " ..
r3 ( score_after_quake - best_score ) )
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 ()
undo.SetUndo ( true )
end
--main ()
xpcall ( main , cleanup )
