Recipe: Quaking Remix 1.4.4

created by Bruno Kestemont

Profile

• Name: Quaking Remix 1.4.4
• ID: 102099
• Shared with: Public
• Parent: None
• Children:
  • Quaking Remix 1.4.0 45 rebuilds
• Created on: November 05, 2016 at 18:01 PM UTC
• Updated on: November 05, 2016 at 18:01 PM UTC
• Description:

  Loop Rebuild modified so that after Remix/Rebuild/Shake/Wiggle etc: it performs some Local Quakes rather than Fuses in an attempt to improve the score.

Code

--Quaking Remix  13/02/2014 Bruno Kestemont based on spvincent
--Outlines: on several zones of X segments, remix, shake, take Y bests, further improve, take best, rebuild it, quake the rebuild result
--Added loop counter and timestamps in output, progress reported in segment Note, changed default ss_control to 2, commented out some print lines
--Added dialog box Bruno Kestemont 02/02/2014
--changed wiggles by recursive wiggles for New Chapter --03/02/2014
--minor adaptations 05/02/2014
--tvdl WORKON dialog added by BK 06/02/2014 
--13/02/2014 fixed dialog box options display
--added remix command, experimental
--19/02/2014 implemented remix on mutable puzzles
--01/03/2014 options skip quake skip rebuild
--26/4/2016 v1.0 shake after each remix
--3/5/2016 v1.1 new remix command used
-- v1.2 fixed remix action !
-- v1.3 loop option on first dialog, only loop as default, user bands management, all remixes considered 16/5/2016
-- v1.3.1 fixed NeverDisableMyBands, disable cut option (it does not work with new remix tool)
-- v1.3.2 fixed recipe name in set note
-- v1.4.0 shakeCI default = 0.31 in wiggleSimple (to be change manually in recipe, not in dialog)
--  and normal shake (CI=1) on a sphere SelectAround(ss,se,radius,nodeselect)
--  rebuild 1 if no remix found (rebuild is a remix of length 3)
--  this version should be slow, but it might gain more
-- added Fast option (skipping all above, only one local shake after each remix)
--v1.4.1 optimized Fast option and adding option for number of best remixes to be refined (low is quick, default 4)
-- optional RebuildOnRemix (option in the code only), if false then rebuild this segments on the credited score
-- but the quake is always on rebuild result (like in Ebola).
--v1.4.2 fixed bug fast
--v1.4.3 jeff101 display, retired fast, adding NoWiggle (wiggle is dangerous for rama map if the AAs are not the right ones)
--v1.4.4 undo.SetUndo
--TO DO ! default all protein (not only loops) & select where to work doen't work !?

recipename="Quaking Remix 1.4.4"

undo.SetUndo(false)

--WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
--OPTIONS FOR DEFAULT SETTINGS to be edited in the code
ss_control =  1
-- ss_control == 1                    Only rebuilds/remix loops
-- ss_control == 2                     Rebuilds everything. Helices and sheets will be temporarily converted to loops
-- ss_control == 3                    Rebuild/remix loops but allows 1+1 non-loop residues in the rebuild segment.
			if ss_control==1 then
				DefaultOption="Default: Only rebuilds loops"
			elseif ss_control==2 then
				DefaultOption="Default: Rebuilds everything converted to loops"
			elseif ss_control==3 then
				DefaultOption="Default: Work on loops with 1 non-loop allowed"
			end
			print(DefaultOption)
			print("Default can easilly be edited in the code")
--of course you can edit other default options bellow
--I could add a dialog box for bands later on, but I suppose default options are ok.
--WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
-- Handy shorts module from tvdl  03/02/2014
normal= (current.GetExplorationMultiplier() == 0)
segCnt=structure.GetCount()
segCnt2=segCnt
while structure.GetSecondaryStructure(segCnt2)=="M" do segCnt2=segCnt2-1 end
-- end Handy shorts module

-- Segment set and list module
-- Notice that most functions assume that the sets are well formed
-- (=ordered and no overlaps)
-- 02-05-2012 TvdL Free to use for non commercial purposes
-- 05-02-2014 BK enhanced (also Free for non commercial purposes)
-- a list is {1,2,3,8,9,10}. A set is {{1,4},{8,10}} or something like that in a 2-dimentional table:
--i	{[1],[2]}
--___________
--1	{ 1 , 4 }
--2	{ 8 , 10}
allowcut= false -- with the new remix tool, I think it's not necessary any more
CUTWORST= false -- not used
shakeCI= 0.31 -- to be changed here manually, not in dialog. Note that there are still shakes at CI=1 occuring
Fast= false -- if true, does only one shake on each remix, and no rebuild if no solution
MaxNbBest= 4  -- we take only the best 4 remixes
NbBest= MaxNbBest -- we take only the best 4 remixes
RebuildOnRemix= true -- rebuild starts on the best remix solution (not on best score)
NoWiggle= false -- if true: never wiggle backbones (protects the rama map), experimental

function cutRange(segFrom, segTo) -- 26/4/2016
	if (allowcut or CUTWORST) and not freeze.IsFrozen(segFrom) then
		if segFrom==segCnt then segFrom=segCnt-1 end
		structure.InsertCut(segFrom)
	end -- else nothing !
	if (allowcut or CUTWORST) and not freeze.IsFrozen(segTo) then
		if segTo==segCnt then segTo=segCnt-1 end
		structure.InsertCut(segTo)
	end -- else nothing !
end

function uncut() -- 18/2/2016
	if (allowcut or CUTWORST) then
		for i=1, segCnt do
			structure.DeleteCut(i)
		end
	end
end

function SegmentExtendSet(set,x,y) -- New BK 05/02/2014 x, y= nb of seg to add before and after each zone
    local result={}
	local x,y=x or 1, y or 1
    for i=1,#set do
		if set[i][1]<1+x then									--first segment
			result[i]={set[i][1],set[i][2]+y}
		elseif set[i][2]>set[#set][2]-y then					--last segment
			result[i]={set[i][1]-x,set[i][2]}
		else
			result[i]={(set[i][1]-x and set[i][1]>x) or set[i][1] , (set[i][2]<=set[#set][2]-y and set[i][2]+y) or set[i][2]}
			result[i]={set[i][1]-x ,set[i][2]+y}
		end
    end
    return result -- overlaps are possible
end

function SegmentListToSet(list) --note: clean duplicates
    local result={}
    local f=0
    local l=-1
    table.sort(list)
    for i=1,#list do
        if list[i] ~= l+1 and list[i] ~= l then
            -- note: duplicates are removed
            if l>0 then result[#result+1]={f,l} end
            f=list[i]
        end
        l=list[i]
    end
    if l>0 then result[#result+1]={f,l} end
    --print("list to set")
    --SegmentPrintSet(result)
    return result
end

function SegmentSetToList(set)
    local result={}
    for i=1,#set do
        --print(set[i][1],set[i][2])
        for k=set[i][1],set[i][2] do
            result[#result+1]=k
        end
    end
    return result
end

function SegmentCleanSet(set) -- clean duplicates
-- Makes it well formed and clean duplicates
    return SegmentListToSet(SegmentSetToList(set))
end

function SegmentInvertSet(set,maxseg)
-- Gives back all segments not in the set
-- maxseg is added for ligand
    local result={}
    if maxseg==nil then maxseg=structure.GetCount() end
    if #set==0 then return {{1,maxseg}} end
    if set[1][1] ~= 1 then result[1]={1,set[1][1]-1} end
    for i=2,#set do
        result[#result+1]={set[i-1][2]+1,set[i][1]-1}
    end
    if set[#set][2] ~= maxseg then result[#result+1]={set[#set][2]+1,maxseg} end
    return result
end

function SegmentInList(s,list)
    table.sort(list)
    for i=1,#list do
        if list[i]==s then return true
        elseif list[i]>s then return false
        end
    end
    return false
end

function SegmentInSet(set,s)
    for i=1,#set do
        if s>=set[i][1] and s<=set[i][2] then return true
        elseif s<set[i][1] then return false
        end
    end
    return false
end

function SegmentJoinList(list1,list2) -- duplicates allowed
    local result=list1
    if result == nil then return list2 end
    for i=1,#list2 do result[#result+1]=list2[i] end
    table.sort(result)
    return result
end

function SegmentJoinSet(set1,set2) -- no duplicates, no overlap
    return SegmentListToSet(SegmentJoinList(SegmentSetToList(set1),SegmentSetToList(set2)))
end

function SegmentCommList(list1,list2) -- what is common to 2 lists
    local result={}
    table.sort(list1)
    table.sort(list2)
    if #list2==0 then return result end
    local j=1
    for i=1,#list1 do
        while list2[j]<list1[i] do
            j=j+1
            if j>#list2 then return result end
        end
        if list1[i]==list2[j] then result[#result+1]=list1[i] end
    end
    return result
end

function SegmentCommSet(set1,set2) -- what is common to set 1 and set 2 (no duplicate)
    return SegmentListToSet(SegmentCommList(SegmentSetToList(set1),SegmentSetToList(set2)))
end

function SegmentSetMinus(set1,set2) -- set1 minus set2 (set1 minus what is common to set 1 and set2)(no duplicate)
    return SegmentCommSet(set1,SegmentInvertSet(set2))
end

function SegmentPrintSet(set)
    print(SegmentSetToString(set))
end

function SegmentSetToString(set)
    local line = ""
    for i=1,#set do
        if i~=1 then line=line..", " end
        line=line..set[i][1].."-"..set[i][2]
    end
    return line
end

function SegmentSetInSet(set,sub)
    if sub==nil then return true end
    -- Checks if sub is a proper subset of set
    for i=1,#sub do
        if not SegmentRangeInSet(set,sub[i]) then return false end
    end
    return true
end

function SegmentRangeInSet(set,range)
    if range==nil or #range==0 then return true end
    local b=range[1]
    local e=range[2]
    for i=1,#set do
        if b>=set[i][1] and b<=set[i][2] then
            return (e<=set[i][2])
        elseif e<=set[i][1] then return false end
    end
    return false
end

function SegmentSetToBool(set)
    local result={}
    for i=1,structure.GetCount() do
        result[i]=SegmentInSet(set,i)
    end
    return result
end
--- End of Segment Set module

-- Module Find Segment Types
function FindMutablesList()
    local result={}
    for i=1,segCnt2 do if structure.IsMutable(i) then result[#result+1]=i end end
    return result
end
function FindMutables()
    return SegmentListToSet(FindMutablesList())
end

function FindFrozenList()
    local result={}
    for i=1,segCnt2 do if freeze.IsFrozen(i) then result[#result+1]=i end end
    return result
end
function FindFrozen()
    return SegmentListToSet(FindFrozenList())
end

function FindLockedList()
    local result={}
    for i=1,segCnt2 do if structure.IsLocked(i) then result[#result+1]=i end end
    return result
end
function FindLocked()
    return SegmentListToSet(FindLockedList())
end

function FindSelectedList()
    local result={}
    for i=1,segCnt do if selection.IsSelected(i) then result[#result+1]=i end end
    return result
end
function FindSelected()
    return SegmentListToSet(FindSelectedList())
end

function FindAAtypeList(aa)
    local result={}
    for i=1,segCnt2 do
        if structure.GetSecondaryStructure(i)== aa then result[#result+1]=i end
    end
    return result
end
function FindAAtype(aa)
    return SegmentListToSet(FindAAtypeList(aa))
end

function FindAminotype(at) --NOTE: only this one gives a list not a set
    local result={}
    for i=1,segCnt2 do
        if structure.GetAminoAcid(i) == at then result[#result+1]=i end
    end
    return result
end
-- end Module Find Segment Types

function SortBest(tab,items) --BACWARD bubble sorting - best on top, only needed items
    for x=1,items do --items do
        for y=x+1,#tab do
            if tab[x][1]<tab[y][1] then
                tab[x],tab[y]=tab[y],tab[x]
            end
        end
    end
    return tab
end

function SelectAround(ss,se,radius,nodeselect)
    if nodeselect~=true then selection.DeselectAll() end
    for i=1, segCnt2 do
        for x=ss,se do
            if structure.GetDistance(x,i)<radius then selection.Select(i) break end
        end
    end
end

--PRIMARY Quaking Rebuild options (renamed)

total_rebuild_loops = 9999 -- exit after rb # exceeds total_rebuild_loops
From_rebuild_length =  9                      -- Maximum rebuild length.
To_rebuild_length =  3                      --  -  Better less than or equal to From_rebuild_length. Suggest not going below 3 for rebuild.
maxFrom_rebuild_length =  11 --28/4/2016
maxTo_rebuild_length =  8--03/02/2014

number_of_rebuilds =  0                   -- default 0 for remix; suggested 8 for the  best score from this number of rebuilds. Change to taste.
										-- if nor remix is found, set at least number_of_rebuilds to 1 if you want a rebuild on this section
number_of_remix=85							--remix is quick here  01/03/2014, we are limited by the 100 available slots
start_idx = 0                                    --  if this is a residue id, use it as the start point. If not
                                                             --  use a  semi-random number computed from the score
--ss_control =  3  			see top for edit recipe default option
                                                                    --  they'll be temporarily converted to loops
loop_forever = true                               --  just keep going. Good for overnight running.
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
kClashImportanceDuringRebuild = 1    -- Might try reducing this to get more a different rebuild than the highest scoring                                                                                                                                        -- one, which might overstate
                                                                            -- the importance of clashing which can be remedied by a quick shake.
-- If, after rebuild, shake, wiggle, etc. , the score is close enough to the starting score, attempt some quakes.
-- close enough is defined by the entry in kQuakingThresholds indexed by the size of segment being rebuilt
kQuakingThresholds = { 8 , 8 , 8 , 15 , 25 , 30 , 35 }
nQuakingThresholds = 7
-- Use of quicksave slots
kOriginalStructureOrNewBest = 1
kBestPostRebuild = 2
--    Stores the  best non-improved structure in Quicksave slot 10 (but can change this below) . The structure that's saved
--    is the best before WS etc, so if this script fails to produce (or even if it does) you can load the structure in
--    slot 10 and try WS in different ways.
kBestNonImprovedStructure = 10
kbestremixslot1=15 -- starting slot, others are 16, 17 ... 100
original_secondary_structure = {}
n_residues = structure.GetCount ()
curr_score = 0
multiplier = 1
-----
-- The Local Quake bit. Copied from LQ1.2 for expediency
-- Bands radiate from a segment to spatially close neighbours
-- Foregoes Fuses and generally tries to optimize for speed

use_ranked_score = 1                       -- if 1, use the score being used for ranking (best normally). If 0, use the
                                                               --  old style scoring system
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                             -- Set to 0 for a single centre. This is the most likely one to change.
opposite_centre = false                       -- overrides second_offset if true
target_wiggle_reduction = 120        -- 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,
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
n_quakes  = 10                                -- A reasonable number 10
band_to_minima = true                -- Band to minima or maxima in the distance. See bands_from_segment for meaning.
do_second_shake = false           -- After creating the bands and pulling with them do a Shake/Wiggle/Shake
                                                          -- if do_second_shake is false (best I think) don't do the last shake.
                                                          -- Speeds things up quite a bit
scale_band_strength = true        -- If true, band strength will be scaled inversely with length
monit = false			-- prints more info

---------- A few more globals
candidate_ids_start = {}
candidate_ids_end = {}
candidate_distances = {}
n_candidates = 0
recent_best_score = 0-- recent best to compare to and see if score was improved
normal= (current.GetExplorationMultiplier() == 0)--03/02/2014

function report_time(start_clock,start_time,clock_msg,time_msg)
  local seconds,minutes,hours,days
  if clock_msg==nil then clock_msg="CPU time" end
  if time_msg==nil then time_msg="Elasped time" end
  print(string.format("%s",os.date()))
  days,remainder=math.modf((os.clock()-start_clock)/(24*60*60))
  hours,remainder=math.modf(remainder*24)
  minutes,remainder=math.modf(remainder*60)
  seconds,remainder=math.modf(remainder*60)
  print(string.format("%s(%02id:%02ih:%02im:%02is)",clock_msg,days,hours,minutes,seconds))
  days,remainder=math.modf(os.difftime(os.time(),start_time)/(24*60*60))
  hours,remainder=math.modf(remainder*24)
  minutes,remainder=math.modf(remainder*60)
  seconds,remainder=math.modf(remainder*60)
  print(string.format("%s(%02id:%02ih:%02im:%02is)",time_msg,days,hours,minutes,seconds))
end

-- 2 Score functions -- NEW from tvd for exploration puzzles BK 03/02/2014
function GetScore(pose)--03/02/2014
    if pose==nil then pose=current end
    local total= pose.GetEnergyScore()
    -- FIX for big negatives

        if total < -999999 and total > -1000001 then total=SegScore(pose) end

    if normal then
        return total
    else
        return total*pose.GetExplorationMultiplier()
    end
end

function SegScore(pose)--03/02/2014
    if pose==nil then pose=current end
    local total=8000
    for i=segStart, segCnt2 do
        total=total+pose.GetSegmentEnergyScore(i)
    end
    return total
end

-- END score functions

global_ci=1 --03/02/2014
function CI(val)--03/02/2014
	global_ci=val
    behavior.SetClashImportance(global_ci)
end

function WiggleSimple(val,how) -- Bruno Kestemont 02/02/2014
	local val, how = val or 1, how or "wa" -- new BK 02/02/2014
		if how == "s" then CI(shakeCI) structure.ShakeSidechainsSelected(1) CI(global_ci)
            elseif how == "wb" then structure.WiggleSelected(val, true, false) -- backbones
            elseif how == "ws" then structure.WiggleSelected(val, false, true) -- sidechains
            elseif how == "wa" then structure.WiggleSelected(val, true, true) -- all
			elseif how== "lw" then structure.LocalWiggleSelected(val) -- local wiggle
        end
end

function WiggleNC(ss, how, iters, minppi)-- 19/8/2016 (ss is useless here), with NoWiggle option
    local valiter=2
    local val=1
    if how==nil then how="wa" end
	if (how=="wa" or how=="lw") and NoWiggle then how="ws" end
    if iters==nil then iters=30 end
	minppi=(recent_best_score-GetScore ())/100 -- 03/02/2014
    if ((minppi==nil) or (minppi<0.001)) then
            minppi=0.001 -- put higher to save time
    end
    if global_ci==1.00 then val=valiter end
    if iters>0 then
		iters=iters-1
        local sp=GetScore()-- 03/02/2014
		WiggleSimple(val,how)
        local ep = GetScore()-- 03/02/2014
        local ig=ep-sp
        if how~="s" then
            if ig > minppi then WiggleNC(ss, how, iters, minppi) end
        end
	how="wa"
    iters=30
    end
end

--[[ 03/02/2014 (old score function)
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 roundto3 ( i )

-- printing convenience

   return i - i % 0.001

end

function roundto6 ( i )

-- printing convenience

   return i - i % 0.000001

end

function  delete_all_bands () -- unless user bands 16/5/2016

    ManageBands ()

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

   most_distant_minimum = segment_ids [ n_segments ]

end

function create_bands ()

  if ( scale_band_strength == true ) then

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

    average_band_length = total_band_length / n_candidates

  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 ()
    --structure.WiggleSelected ( 1 )
	WiggleNC(seg, "wa", 30) -- 02/02/2014

    score_after_wiggle = GetScore ()

    score_reduction_after_wiggle = init_score - score_after_wiggle

    -- Adaptively change band strength based on how much the previous wiggle overshot/undershot the target. Could
    -- be made much more elaborate.

    if ( score_reduction_after_wiggle > target_wiggle_reduction * 1.2 ) then
        band_strength = band_strength * ( 0.85 + 0.1 *pseudorandom ( score_after_wiggle ) )
     elseif ( score_reduction_after_wiggle  < target_wiggle_reduction * 0.8 ) then
        band_strength =  band_strength + ( 1 - band_strength ) * ( 0.05 + 0.1 *pseudorandom ( score_after_wiggle ) )
    end

    delete_all_bands ()

    structure.ShakeSidechainsSelected ( 1 )

    gas = GetScore () - score_after_wiggle

    if ( monit ) then
       print ( "law " .. roundto3 ( init_score - score_after_wiggle ) .. " gas " ..  roundto3 ( gas ) .. " band_strength " .. roundto3 ( band_strength )  )
    end

    --structure.WiggleSelected ( 8 )
	WiggleNC(seg, "wa", 30) -- 02/02/2014

    if ( do_second_shake == true ) then

      score_a2w = GetScore ()
      structure.ShakeSidechainsSelected ( 1 )
      score_a2s = GetScore ()

      if ( score_a2s - score_a2w  > 1 ) then
          --structure.WiggleSelected  ( 4 )
		  WiggleNC(seg, "wa", 30) -- 02/02/2014
      end
   end

end

function local_quake ( start_idx )

  band_strength = 0.5
  CI ( 1 )

  recentbest.Save  ()

  best_score_lq = GetScore ()

  inc = Coprime ( n_residues )

  idx_lq = start_idx
  print("Quaking from seg "..start_idx.." and score: "..roundto3(curr_score))

  for i =  1, n_quakes  do

    recentbest.Restore  ()

     -- Need to record any improvements from the previous iteration,.

    score_lq = GetScore()
    if ( score_lq > best_score_lq ) then
       best_score_lq = score_lq
--     print ( "      Gain to " .. roundto3 ( score_lq - score_lq % 0.01 ) )
    end

    delete_all_bands ()

    n_candidates = 0 -- ??
    get_band_segment_data  ( idx_lq )

    if ( opposite_centre == true ) then
      second_idx = most_distant_minimum
      get_band_segment_data  ( most_distant_minimum )
     elseif ( second_offset ~= 0 ) then

      second_idx = idx_lq + 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 -- ?? it's always zero in this loop ??!?

      create_bands ()

      if ( second_offset ~= 0 ) then
--       print ( "    i: " .. i .. " seg ".. idx_lq .. "," .. second_idx )
       else
--       print ( "    i: " .. i .. " seg ".. idx_lq )
      end

     Quake ( idx_lq )

   end

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

  end -- for

  recentbest.Restore  ()

  delete_all_bands ()

end
------
function IsSegmentALoop ( first , last )

  -- Returns number of non-loop  residues (helices or sheets)

   n_non_loops = 0
   for i = first , last do
       if ( structure.GetSecondaryStructure ( i ) ~= "L" ) then
          n_non_loops = n_non_loops + 1
        end
     end

   return  n_non_loops

end

function ConvertToLoop ( first , last )

    for  k = first , last do
        if ( original_secondary_structure [ k ] ~= "L" ) then
            selection.DeselectAll ( )
            selection.Select ( k )
            structure.SetSecondaryStructureSelected ( "L" )
        end
    end
end

function RestoreStructure ( first , last )

    for k = first, last do
        selection.DeselectAll ()
        selection.Select ( k )
        structure.SetSecondaryStructureSelected  ( original_secondary_structure [ k ] )
    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 nudge (  score )

    selection.SelectAll ()
    --recentbest.Save ()
    CI ( 0.2  + 0.2 * pseudorandom ( score ) )
    --structure.WiggleSelected ( 1 )
	WiggleNC(seg, "wa", 2) -- 02/02/2014
    CI ( 1 )
   -- structure.WiggleSelected ( 8 )
	WiggleNC(seg, "wa", 30) -- 02/02/2014
    recentbest.Restore ()
end

function  WiggleShakeWiggle ( start_idx , end_idx )

     -- returns the score after doing a WSW:  shortcircuiting if things don't look promising

    score_a = GetScore ()

    --structure.WiggleSelected ( 15 )
	WiggleNC(seg, "wa", 30) -- 02/02/2014
    recentbest.Restore ()

    score_after_wiggle = GetScore ()

    if ( score_after_wiggle - score_a  < 10*multiplier and recent_best_score - score_after_wiggle > 30*multiplier  ) then
        --    The wiggle got stuck and didn't achieve anything
        nudge ( score_after_wiggle )
        score_after_wiggle = GetScore ()
    end

    if ( monit ) then
        --print (  "After first wiggle (bb) : "..  roundto3 ( recent_best_score - score_after_wiggle ) )
    end

    if ( recent_best_score - score_after_wiggle  > 100*multiplier ) then
        -- Not worth continuing
       return  score_after_wiggle
    end

    structure.ShakeSidechainsSelected ( 1 )

    score_after_second_shake = GetScore ()

    if ( true ) then --score_after_second_shake - score_after_wiggle > 0.5 ) then
        --structure.WiggleSelected ( 10 )
		WiggleNC(seg, "wa", 30) -- 02/02/2014
        recentbest.Restore ()
        score_after_second_wiggle = GetScore  ()

        if ( true ) then -- score_after_second_wiggle - score_after_second_shake  < 1  and recent_best_score  - score_after_second_wiggle < 50*multiplier ) then

            -- The second shake did something but the subsequent wiggle didn't

            nudge ( score_after_second_wiggle )

            score = GetScore  ()

            return  score
        else
            return  score_after_second_wiggle
        end
    else
        return score_after_second_shake
    end

end

function  post_rebuild (  score_after_rebuild , start_idx , end_idx )

    --  Tries to recover after a rebuild.

    if ( recent_best_score - score_after_rebuild > 3000*multiplier  ) then
       structure.ShakeSidechainsSelected (  2 )
    elseif  ( recent_best_score - score_after_rebuild > 1500*multiplier  ) then
        structure.ShakeSidechainsSelected ( 1 )
    end

    --    First try  WSW

    curr_score = WiggleShakeWiggle ( start_idx , end_idx )

end

function updatebest(curr_score)-- 22/02/2014 and print improvements
    if ( curr_score > recent_best_score ) then
		local gain=curr_score - recent_best_score
		print ( "Improvement to "..  roundto3(curr_score) .. ". This gain "..roundto6(gain))
		structure.SetNote(note_number,string.format("(%s) %.3f + %s (%i) %.3f",user.GetPlayerName(), starting_score, recipename, rb,roundto3(curr_score)))
		jeff101Display(recent_best_score, curr_score, rb)
		
		recent_best_score = curr_score -- note that we work only with "recent best" - TO DO: verify

		improvement_made = true

		save.Quicksave(kOriginalStructureOrNewBest)
 
    end
end

function updatebestRemix(curr_score) -- not to be used too early (kbestremixslot1 is needed on start remix)
	if curr_score>bestremixscore then
		bestremixscore=curr_score
		save.Quicksave(kbestremixslot1) -- we keep here the very best remix
	end
end					
					
function checkLoop(start_idx ,  end_idx)
	if IsSegmentALoop ( start_idx , end_idx ) > 0 then return false end
	return true
end

function AfterRemix( start_idx  ,  end_idx ) -- on selected
		recentbest.Save()
		WiggleNC(ss, "ws", 2, 0.1)
		WiggleNC(ss, "s", 1, 0.1) -- with CI= 0.31
		WiggleNC(ss, "lw", 1, 0.1) -- it's very slow !
		--WiggleNC(ss, "wa", 2, 0.1) -- selected
		selection.DeselectAll ()
		SelectAround(start_idx,end_idx,9)
		structure.ShakeSidechainsSelected(1) -- on a sphere with CI=1
		recentbest.Restore()
		selection.DeselectAll ()
		selection.SelectRange (  start_idx  ,  end_idx  )
end
				
function rebuild_seg ( start_idx , end_idx ) -- remix is here (with or without additional rebuild)
	local rebuild_lengthseg = end_idx - start_idx + 1 --22/02/2014 initialization 
	local remixscores= {}
	local remixdone= false
--      return codes
--
--      0       rebuild had no effect: seg probably locked
--      1       no improvement
--      2       improvement

  selection.DeselectAll ()
  selection.SelectRange (  start_idx  ,  end_idx  )

  -- for some reason, do_local_rebuild ( 1)  often leaves  the protein unchanged
  -- Compute a value delta which will be used later to determine if we are
  -- "sufficiently" close to our original  best score.  If it turns out later that
  -- we are, then nudge it around

  rebuild_length = end_idx - start_idx + 1 -- duplicate but it's safe because we need it here (could be local)

  if ( rebuild_length > nQuakingThresholds ) then
    delta = kQuakingThresholds [ nQuakingThresholds ] * multiplier       -- sanity check for lengthy rebuild segs
  else
    delta = kQuakingThresholds [ rebuild_length ] * multiplier -- multiplier = 1 most of the time
  end
  
  
  k  = 0
  repeat
     k = k + 1
	 if do_remix and rebuild_lengthseg>2 and  rebuild_lengthseg< 10 and ss_control~=2 then -- 16/5/2016
		save.LoadSecondaryStructure () --19/02/2014
		print("Remixing",start_idx,"-",end_idx,". Length=",rebuild_lengthseg)
		recentbest.Save () --19/02/2014
		cutRange( start_idx -1  ,  end_idx +1) --26/4/2016
		save.Quicksave(kbestremixslot1) -- starting and best remix score will be here
		--save.Quicksave(kstartremixslot)
		
		bestremixscore=GetScore ()-900000 -- threshold to be adapted
		remixdone= false
		remixscores= {}
		local nfound=structure.RemixSelected(kbestremixslot1, number_of_remix) -- saves different remix results in max x (10?) slots
		print("Nr of remix found "..nfound)
		if nfound ~= 0 then --16/5/2016
			if nfound >= MaxNbBest then NbBest= MaxNbBest else NbBest=nfound end
			noRemix=false
			for i=1 , nfound do --15/5/2016
				save.Quickload(kbestremixslot1+i-1)
				--structure.RemixSelected()

				structure.ShakeSidechainsSelected(1)

				curr_score = GetScore ()
				remixscores[#remixscores+1]={curr_score, i}
				--print("Score after remix ",i,":",roundto3(curr_score))
				save.Quicksave(kbestremixslot1+i-1) -- replace the old value by this shaked one
			end
				SortBest(remixscores,NbBest) -- sort only to find the needed best
				bestremix=0
				print ("Selecting", NbBest, "best remixes:")
				for i= 1, NbBest do
					save.Quickload(kbestremixslot1+remixscores[i][2]-1)
					print("Remix ", remixscores[i][2], "score:", roundto3(remixscores[i][1]))
					if not Fast then
						AfterRemix( start_idx  ,  end_idx )
					end
					curr_score = GetScore ()
					if curr_score>bestremixscore then
						bestremixscore=curr_score 
						bestremix= remixscores[i][2]
						save.Quicksave(kbestremixslot1) -- we keep here the very best remix
					end
				end
			save.Quickload(kbestremixslot1) -- We take ONLY the best one (TO DO: second best etc ?)
			uncut()--26/4/2016
			curr_score= GetScore ()
			if curr_score ~= recent_best_score then --16/5/2016
				updatebestRemix(curr_score)
				structure.ShakeSidechainsSelected(1)
				curr_score = GetScore ()
				print("Best score after remix round and fix: "..roundto3(curr_score).." on remix "..bestremix)
				updatebestRemix(curr_score)
				save.Quickload(kbestremixslot1) 
				WiggleNC(ss, "wa", 2, 0.1)
				curr_score = GetScore ()
				updatebestRemix(curr_score)
				save.Quickload(kbestremixslot1) 
				updatebest(curr_score) -- 22/02/2014 and print improvements
			end
		else --16/5/2016
			uncut()
			save.Quickload(kOriginalStructureOrNewBest)
			noRemix=true
		end
	 end
     score_after_remix = GetScore ()
   until  score_after_remix ~= recent_best_score or k  > 2

   if ( score_after_remix == recent_best_score ) then

      -- If  the scores are  still equal  we probably have a locked segment
	print("Remix didn't change anything")
     return 0 -- WARNING then no rebuild, no quake !?! TO DO: evaluate or optional ?

  end
  if not RebuildOnRemix or score_after_remix <  recent_best_score - delta    then
	  recentbest.Restore () -- 22/02/2014 further rebuilds not on loosing remixed
  end
  	  recentbest.Save ()
	  score_before_rebuild = GetScore ()
  if (  number_of_rebuilds > 0 ) then
	print("Rebuilding "..start_idx.."-"..end_idx.." on score "..roundto3(score_before_rebuild))
    CI ( kClashImportanceDuringRebuild  )
	ConvertToLoop ( start_idx ,  end_idx )--
    structure.RebuildSelected ( number_of_rebuilds)
    CI ( 1 )
    recentbest.Restore ()
	structure.ShakeSidechainsSelected(1) -- with CI=1 on the best rebuild only
	recentbest.Restore ()
  end

  -- Now do the post rebuild wiggle and shakes

  score_after_rebuild = GetScore ()

  selection.SelectAll ()



   save.Quicksave ( kBestPostRebuild )

    post_rebuild ( score_after_rebuild , start_idx , end_idx  ) -- it calculates also the resulting curr_score

    if (  recent_best_score - curr_score < delta ) then
      print ( "   cur score " .. roundto3 ( curr_score ) )

        local_quake ( start_idx )
    end

    --if ( ss_control > 1 ) then
      save.LoadSecondaryStructure ()
    --end

    curr_score = GetScore ()

	updatebest(curr_score) -- 22/02/2014 and print improvements


  if (  curr_score > qs10_best_score  and improvement_made == false ) then

   recentbest.Save ()
   qs10_best_score = curr_score
   save.Quickload ( kBestPostRebuild )

    -- Bracket the area with frozen segments for easy recognition
    selection.DeselectAll ()

    if ( start_idx > 1 ) then
      selection.Select  ( start_idx - 1 )
    end

    if ( end_idx <  n_residues ) then
      selection.Select  ( end_idx + 1 )
    end

    freeze.FreezeSelected ( true , false )

    save.Quicksave ( kBestNonImprovedStructure )
    recentbest.Restore ()

  end

  if ( improvement_made == true ) then
    return 2
  else
    return 1
  end

end

function rebuild ( rebuild_length  )
    n_possible_segs = n_residues - rebuild_length + 1 -- max number of potential zones
    inc = Coprime ( n_possible_segs  ) -- mostly 1, or 31 or 37 don't ask me why
    --    The point of  this Coprime, start_idx business is so we go through the protein in  a
    --    nonuniform way. Rather than rebuilding 1-4, 2-5 , 3-6 for example, we might do 21 4, 53-57, 3-6
    --    or something like  that
	-- example: 100 residues, length 9 => n_possible_segs= 91
	-- inc = 31
	-- it tries 91 times from seg 1, 32, 64, 112-91= 20, 51, 82, 113-91= 22, 53, 84, 115-91= 24 etc
    for  i = 1 , n_possible_segs do -- it tries everything but it verifies afterwards if it is in the right zone (WORKONbool)
        improvement_made = false -- global updated in rebuild_seg (similar to k)
        end_idx = start_idx + rebuild_length - 1

		if WORKONbool[start_idx] and WORKONbool[end_idx]  then -- 05/02/2014
			if (  start_idx >= min_residue and  end_idx <= max_residue and
						( ss_control == 2 or
					   ( ss_control == 1 and IsSegmentALoop ( start_idx , end_idx  ) == 0 ) or -- nb loops =0
					   ( ss_control == 3 and IsSegmentALoop ( start_idx , end_idx  ) < 2  ) ) ) then
				rb = rb + 1
				if rb > total_rebuild_loops then return end
				print (string.format("rb %i/%i %i-%i (%i/%i) Score: %.3f Best score: %.3f",rb,total_rebuild_loops,start_idx,end_idx,i,n_possible_segs,current.GetScore(),recent_best_score))
				if rb%10==0 then report_time(start_clock,start_time) end
				save.Quickload ( kOriginalStructureOrNewBest )

				if ( ss_control > 1 ) then
					ConvertToLoop ( start_idx ,  end_idx )
				end
				recent_best_score = GetScore () -- reinitialises the recent best here

					k = rebuild_seg ( start_idx, end_idx ) -- k= 1 if no improvement, 2 if improvement, 0 if rebuild failure
					print("0=failure; 2=improvement:",k)

			  --if ( ss_control > 1 ) then
				   save.LoadSecondaryStructure ()
				--end
			end
		end
      start_idx = start_idx + inc -- look for the next possible zone
      if ( start_idx >  n_possible_segs ) then
        repeat
          start_idx =  start_idx - n_possible_segs
        until start_idx<=n_possible_segs -- starts again from the first segments
      end
    end        --    for i
end


start_clock=os.clock()
start_time=os.time()
print(string.format("%s",os.date()))

note_number=structure.GetCount()
for seg=structure.GetCount(),1,-1 do
  if structure.GetNote(seg)~="" then break end
  note_number=seg
end
print(string.format("Recording results in Note for segment %i",note_number))
starting_score=current.GetScore()

--archiving user bands
ubandlist={}
USERBANDS=false
function archivebands() -- to keep user bands
	--if you want to keep your own bands, for example between well formed sheets
	ubcount=band.GetCount()	-- global
	
	for i=1, ubcount do
		ubandlist[i]=i
	end
	if #ubandlist==0 then
		USERBANDS=false
	else
		USERBANDS=true
	end
	return ubandlist
	
end

NeverDisableMyBands= true

function ManageBands() --if bands can't be saved, prevent qr from deleting them
	if ubcount==0 or ubcount==nil then band.DeleteAll()
	else
		local bands=band.GetCount()
		if bands>ubcount then
			for i=bands, ubcount+1, -1 do band.Delete(i) end -- TO DO: il reste peut etre une bande ici
		end
	end
	if not NeverDisableMyBands then band.DisableAll() end
end

-- Tidy up if necessary
archivebands()

ManageBands ()
n_residues = structure.GetCount ()

--if ( ss_control  > 1 ) then
    for i = 1, n_residues do
        original_secondary_structure [ i ] = structure.GetSecondaryStructure ( i )
    end
--end

save.SaveSecondaryStructure ()
rb = 0
save.Quicksave ( kOriginalStructureOrNewBest )
CI ( 1 )
recent_best_score = GetScore ()
print (  "Start score : " .. roundto3 ( recent_best_score ) )

if ( recent_best_score  == 0 ) then
    print ( "Suggest quitting as starting score < = 0" )
end

function print_ss_control()
	if ( ss_control == 1 ) then
		print (  "Loops only" )
	  elseif ( ss_control == 2 ) then
		print (  "Loops, sheets, and helices" )
	 else
	   print ( "Loops with one non-loop allowed" )
	end
end

if ( max_residue > n_residues ) then
    max_residue = n_residues
end

if ( min_residue < 1 ) then
    min_residue = 1
end

if ( max_residue <= min_residue ) then
    print ( "Rebuild  range error" )
    loop_forever = false
end

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

if ( true ) then -- min_residue  > 1 or  max_residue < n_residues ) then
    print  ( "Rebuild  range " .. min_residue .. " to " .. max_residue )
   else
    print ( "Rebuild all" )
end

qs10_best_score = 0

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

--START jeff101 display (copy this after the starting score and edit it)
function round(x)--cut all afer 3-rd place
    return x-x%0.001
end
function down(x)
    return x-x%1
end
ssstr=round(recent_best_score) -- init jeff101 function 19/8/2016
function jeff101Display(genScore, genEndScore, gen)
        local genGain=genEndScore-genScore
        if gen>1 and (gen-1)==5*down((gen-1)/5) then
            ssstr=(ssstr..'\n  ') -- start new line every 5 generations
        end
        if genGain>=0 then
            ssstr=(ssstr..' + '..round(genGain))
        else
            ssstr=(ssstr..' - '..(-round(genGain)))
        end
        print(ssstr..' = '..round(genEndScore))
end
--END jeff101 display

function AskOptions()--03/02/2014
	WORKON={{1,structure.GetCount()}} -- all segments 06/02/2014
	if ss_control~=2 then								--06/02/2014
		WORKON=SegmentSetMinus(WORKON,FindAAtype("E"))
		WORKON=SegmentSetMinus(WORKON,FindAAtype("H"))
	end
	if ss_control==3 then
		WORKON=SegmentExtendSet(WORKON,1,1)
	end
	local ask=dialog.CreateDialog("Quaking remix Options")
	repeat
		ask.l1=dialog.AddLabel("Remix_length:")
		ask.From_rebuild_length=dialog.AddSlider("From", From_rebuild_length, 1, maxFrom_rebuild_length, 0)
		ask.To_rebuild_length=dialog.AddSlider("To", To_rebuild_length, 1, maxTo_rebuild_length, 0)
		ask.number_of_remix=dialog.AddSlider("Max nb of remix",number_of_remix, 0, 85, 0)
		--ask.allowcut=dialog.AddCheckbox("Cut before remix ",allowcut)
		--ask.Fast=dialog.AddCheckbox("Fast (only 1 local shake after remix)",Fast)
		ask.MaxNbBest=dialog.AddSlider("Work on best:",MaxNbBest, 0, number_of_remix, 0)
		ask.NoWiggle=dialog.AddCheckbox("Never wiggle backbones",NoWiggle)
		ask.number_of_rebuilds=dialog.AddSlider("Nb of rebuilds",number_of_rebuilds, 0, 30, 0)
		ask.n_quakes=dialog.AddSlider("Nb of quakes",n_quakes, 0, 15, 0)

		ask.ss_control=dialog.AddCheckbox("Include 1 seg before & after loops",false) -- 16/5/2016
		ask.SEL= dialog.AddCheckbox("(Re)select where to work on ",false) -- 05/02/2014
		ask.second_offset=dialog.AddSlider("Max Nb centres",second_offset, 0, 3, 0)--not sure it's useful
		if USERBANDS then
			ask.NeverDisableMyBands= dialog.AddCheckbox("Never disable my bands",true) -- 16/5/2016
		end
		ask.ok = dialog.AddButton("OK", 1)
		ask.cancel = dialog.AddButton("Cancel", 0)

		askresult=dialog.Show(ask)
		if askresult > 0 then
        To_rebuild_length=ask.To_rebuild_length.value
        From_rebuild_length=ask.From_rebuild_length.value
		second_offset=ask.second_offset.value
		n_quakes=ask.n_quakes.value
		number_of_rebuilds=ask.number_of_rebuilds.value
		--Fast=ask.Fast.value
		NoWiggle=ask.NoWiggle.value
		--if HASMUTABLE then--19/02/2014
			--do_remix=ask.remix.value --13/02/2014
		number_of_remix=ask.number_of_remix.value
		MaxNbBest=ask.MaxNbBest.value
		if number_of_remix> 0 then do_remix=true else do_remix=false end
		--allowcut=ask.allowcut.value
		if USERBANDS then
			NeverDisableMyBands=ask.NeverDisableMyBands.value -- 16/5/2016
		end
		if ask.ss_control.value then ss_control=3 end
		--end
		if ask.SEL.value then 											-- 05/02/2014
            local SelMode={}
            SelMode.askignorefrozen=false
            SelMode.defignorefrozen=true
            SelMode.askignorelocks=false
            SelMode.defignorelocks=true
            SelMode.askligands=false
            SelMode.defligands=false
            WORKON=AskForSelections("tvdl Select where to work on",SelMode) -- Reinitialization of WORKON
            print("Selection is now, reselect if not ok:")
            print(SegmentSetToString(WORKON))							-- 05/02/2014
			ask.l4=dialog.AddLabel("Selection is now, reselect if not ok:")
			ask.l5=dialog.AddLabel(SegmentSetToString(WORKON))
            if askresult==1 then askresult=4 end --to force return to main menu
        end
		        -- Do not try to rebuild frozen or locked parts or ligands
        WORKON=SegmentSetMinus(WORKON,FindFrozen()) -- duplicates will be removed here
        WORKON=SegmentSetMinus(WORKON,FindLocked())
        WORKON=SegmentSetMinus(WORKON,FindAAtype("M"))
		end
	until askresult<2
	return askresult > 0
end

-- Module AskSelections
-- 02-05-2012 Timo van der Laan, Free to use for non commercial purposes

function AskForSelections(title,mode)									-- 05/02/2014
    local result={{1,structure.GetCount()}} -- All segments, reinitialize always
    if mode == nil then mode={} end 							--what will be in the dialog or not
    if mode.askloops==nil then mode.askloops=true end
    if mode.asksheets==nil then mode.asksheets=true end
    if mode.askhelixes==nil then mode.askhelixes=true end
	if mode.askExtendZones==nil then mode.askExtendZones=true end	-- 05/02/20
    if mode.askligands==nil then mode.askligands=false end
    if mode.askselected==nil then mode.askselected=false end
    if mode.asknonselected==nil then mode.asknonselected=false end
    if mode.askmutateonly==nil then mode.askmutateonly=false end
    if mode.askignorelocks==nil then mode.askignorelocks=false end
    if mode.askignorefrozen==nil then mode.askignorefrozen=false end
    if mode.askranges==nil then mode.askranges=true end
    if mode.defloops==nil then mode.defloops=true end -- default settings on screen
    if mode.defsheets==nil then mode.defsheets=false end
    if mode.defhelixes==nil then mode.defhelixes=false end
	if mode.defExtendZones==nil then mode.defExtendZones=true end	-- 05/02/20
    if mode.defligands==nil then mode.defligands=false end
    if mode.defselected==nil then mode.defselected=false end
    if mode.defnonselected==nil then mode.defnonselected=false end
    if mode.defmutateonly==nil then mode.defmutateonly=false end
    if mode.defignorelocks==nil then mode.defignorelocks=false end
    if mode.defignorefrozen==nil then mode.defignorefrozen=false end
	if ss_control==1 then -- note: default ss_control must be initialized in the code above
		mode.defloops, mode.defExtendZones=true, false
	elseif ss_control==2 then
		mode.defloops, mode.defsheets, mode.defhelixes, mode.defExtendZones=true, true, true, false
	elseif ss_control==3 then
		mode.defloops, mode.defExtendZones=true, true
	end
    local Errfound=false
  repeat
    local ask = dialog.CreateDialog(title)
    if Errfound then
        ask.E1=dialog.AddLabel("Try again, ERRORS found, check output box")
        result={{1,structure.GetCount()}} --reset start
        Errfound=false
    end
    if mode.askloops then
        ask.loops = dialog.AddCheckbox("Work on loops",mode.defloops)
    elseif not mode.defloops then
        ask.noloops= dialog.AddLabel("Loops will be auto excluded")
    end
    if mode.askhelixes then
        ask.helixes = dialog.AddCheckbox("Work on helixes",mode.defhelixes)
    elseif not mode.defhelixes then
        ask.nohelixes= dialog.AddLabel("Helixes will be auto excluded")
    end
    if mode.asksheets then
        ask.sheets = dialog.AddCheckbox("Work on sheets",mode.defsheets)
    elseif not mode.defsheets then
        ask.nosheets= dialog.AddLabel("Sheets will be auto excluded")
    end
    if mode.askExtendZones then
        ask.ExtendZones = dialog.AddCheckbox("Include 1 seg before & after structure",mode.defExtendZones)
    elseif not mode.defExtendZones then
        ask.noExtendZones= dialog.AddLabel("No additional seg")
    end	
    if mode.askligands then
        ask.ligands = dialog.AddCheckbox("Work on ligands",mode.defligands)
    elseif not mode.defligands then
        ask.noligands= dialog.AddLabel("Ligands will be auto excluded")
    end
    if mode.askselected then ask.selected = dialog.AddCheckbox("Work only on selected",mode.defselected) end
    if mode.asknonselected then ask.nonselected = dialog.AddCheckbox("Work only on nonselected",mode.defnonselected) end
    if mode.askmutateonly then ask.mutateonly = dialog.AddCheckbox("Work only on mutateonly",mode.defmutateonly) end
    if mode.askignorelocks then
        ask.ignorelocks =dialog.AddCheckbox("Dont work on locked ones",true)
    elseif mode.defignorelocks then
        ask.nolocks=dialog.AddLabel("Locked ones will be auto excluded")
    end
    if mode.askignorefrozen then
        ask.ignorefrozen = dialog.AddCheckbox("Dont work on frozen",true)
    elseif mode.defignorefrozen then
        ask.nofrozen=dialog.AddLabel("Frozen ones will be auto excluded")
    end
    if mode.askranges then
        ask.R1=dialog.AddLabel("Or put in segmentranges. Above selections also count")
        ask.ranges=dialog.AddTextbox("Ranges","")
    end
    ask.OK = dialog.AddButton("OK",1) ask.Cancel = dialog.AddButton("Cancel",0)
    if dialog.Show(ask) > 0 then
        -- We start with all the segments including ligands
        if mode.askloops then mode.defloops=ask.loops.value end
        if not mode.defloops then
            result=SegmentSetMinus(result,FindAAtype("L"))
        end
        if mode.asksheets then mode.defsheets=ask.sheets.value end
        if not mode.defsheets then
            result=SegmentSetMinus(result,FindAAtype("E"))
        end
        if mode.askhelixes then mode.defhelixes=ask.helixes.value end
        if not mode.defhelixes then
            result=SegmentSetMinus(result,FindAAtype("H"))
        end
	    if mode.askExtendZones then mode.defExtendZones=ask.ExtendZones.value end -- BK 05/02/2014
        if mode.defExtendZones then
            result=SegmentExtendSet(result,1,1) -- will extend loops etc with 1 seg before and 1 after
        end
        if mode.askligands then mode.defligands=ask.ligands.value end
        if not mode.defligands then
            result=SegmentSetMinus(result,FindAAtype("M"))
        end
        if mode.askignorelocks then mode.defignorelocks=ask.ignorelocks.value end
        if mode.defignorelocks then
            result=SegmentSetMinus(result,FindLocked())
        end
        if mode.askignorefrozen then mode.defignorefrozen=ask.ignorefrozen.value end
        if mode.defignorefrozen then
            result=SegmentSetMinus(result,FindFrozen())
        end
        if mode.askselected then mode.defselected=ask.selected.value end
        if mode.defselected then
            result=SegmentCommSet(result,FindSelected())
        end
        if mode.asknonselected then mode.defnonselected=ask.nonselected.value end
        if mode.defnonselected then
            result=SegmentCommSet(result,SegmentInvertSet(FindSelected()))
        end
        if mode.askranges and ask.ranges.value ~= "" then
            local rangetext=ask.ranges.value
            local function Checknums(nums)
                -- Now checking
                if #nums%2 ~= 0 then
                    print("Not an even number of segments found")
                    return false
                end
                for i=1,#nums do
                    if nums[i]==0 or nums[i]>structure.GetCount() then
                        print("Number "..nums[i].." is not a segment")
                        return false
                    end
                end
                return true
            end

            local function ReadSegmentSet(data)
                local nums = {}
                local NoNegatives='%d+' -- - is not part of a number
                local result={}
                for v in string.gfind(data,NoNegatives) do
                    table.insert(nums, tonumber(v))
                end
                if Checknums(nums) then
                    for i=1,#nums/2 do
                        result[i]={nums[2*i-1],nums[2*i]}
                    end
                    result=SegmentCleanSet(result)
                else Errfound=true result={} end
                return result
            end
            local rangelist=ReadSegmentSet(rangetext)
            if not Errfound then
                result=SegmentCommSet(result,rangelist)
            end
        end
    end
  until not Errfound
    return result
end
-- end of module AskSelections
WORKON={{1,segCnt2}}  													-- for To Do
WORKONbool={}
function InitWORKONbool()
    WORKONbool=SegmentSetToBool(WORKON)
end

function ZoneTest() --05/02/2014
	maxlenght=0--the longest zone
	local currenlenght=0
	for i= min_residue, max_residue do
	    if WORKONbool[i] then
			currenlenght = currenlenght + 1
			if currenlenght>maxlenght then
				maxlenght=currenlenght
			end
		else
			currenlenght=0--end of this zone, reset for further zone
		end
	end
	print("Longer zone=",maxlenght)
	return maxlenght
end

function DumpErr(err) --03/02/2014
	start,stop,line,msg=err:find(":(%d+):%s()")
	err=err:sub(msg,#err)
	p('---')
	if err:find('Cancelled')~=nil then
		p("User stop.")
	else
		p("unexpected error detected.")
		p("Error line:", line)
		p("Error:", err)
	end
	LastWish()
end

function LastWish()--03/02/2014
	undo.SetUndo(true)
	print("Restoring CI, best result and structures")
	CI(1)
	save.Quickload ( kOriginalStructureOrNewBest )
	recentbest.Restore()
end

MutList=FindMutablesList()--19/02/2014 NOT USED
HASMUTABLE= (#MutList>0)
if HASMUTABLE then print("Mutables found") end
	
function MAIN()

	InitWORKONbool() -- WORKON is known from call, here list of available segments
	print_ss_control()
	--if ss_control==(1 or 3) then
		local check_length=ZoneTest()
		if check_length < To_rebuild_length or check_length < From_rebuild_length then
			print("Adjusting to available length")
			if check_length < To_rebuild_length then
				To_rebuild_length=check_length
				maxTo_rebuild_length = check_length
			end
			if check_length < From_rebuild_length then
				From_rebuild_length=check_length
				maxFrom_rebuild_length = check_length
			end
		end
	--end

		if From_rebuild_length<To_rebuild_length then increment=1 else increment=-1 end
		repeat
			for rebuild_length = From_rebuild_length , To_rebuild_length , increment  do
				 rebuild ( rebuild_length )
				 if rb > total_rebuild_loops then break end
			end
			if rb > total_rebuild_loops then break end
		until  loop_forever == false
end

if AskOptions() then
    xpcall(MAIN,DumpErr)
end

save.Quickload ( kOriginalStructureOrNewBest )

report_time(start_clock,start_time)

Comments

[Bruno Kestemont]

Optimized:
Remix & shake
Takes the 4 best ones: local shake, wiggle sidechains, local wiggle selected
Take this best one and further fix with shake, wiggle

[Bruno Kestemont]

–v1.4.1 adding option for number of best remixes to be refined (low is quick, default 4)
– optional RebuildOnRemix (option in the code only), if false then rebuild this segments on the credited score
– but the quake is always on rebuild result (like in Ebola).
–v1.4.3 jeff101 display, retired fast, adding NoWiggle (wiggle is dangerous for rama map if the AAs are not the right ones)
–v1.4.4 undo.SetUndo for better performance

[Bruno Kestemont]

Default all protein (not only loop), better for "no rebuild" puzzles. Change is if you want to idealize loops.

[Mike Cassidy]

I am remixing 3-9 segs in that order even though I set it manually from 9-4 segs.

[Bruno Kestemont]

I'll correct it asap

[Bruno Kestemont]

But if existing loops are smaller, the receipe will adapt.