Recipe: General Assembly v1.0

created by Serca

Profile

• Name: General Assembly v1.0
• ID: 108947
• Shared with: Public
• Parent: None
• Children: None
• Created on: December 21, 2024 at 09:44 AM UTC
• Updated on: December 21, 2024 at 11:24 AM UTC
• Description:

  This tiny script is designed to Rebuild the entire protein in a complex.

  It is specifically intended for rebuilds with no fusing at all or only very brief fusing. The goal is to minimize protein rigidity and avoid mojo by keeping the use of WiggleAll to a minimum.

  The script performs well in the early puzzle stages, making the backbone more natural and stable. It helps lock in the protein structure right after manual folding, before switching to standard DRW recipes. ; Its approach is to generate a large number of solutions (hundreds or even thousands) and evaluate them using the original scoring system without performing a Fuse.

  There are two key settings: ; 1) “Solutions to search” ; The number of Rebuild solutions generated in a single run (one pass of the protein). ; The default value is twice the number of protein segments, i.e., 100–300 rebuilds. ; 2) “Make Fuse” (slider): ; Controls how often fusing happens. Specifically, it triggers a Fuse after every X iterations that don’t score points. ; This is the easiest way to adjust the depth and speed of the Rebuild.

Code

--[[ General Assembly

This tiny script is designed to Rebuild the entire protein in a complex.

It is specifically intended for rebuilds with no fusing at all or only very brief fusing. The goal is to minimize protein rigidity and avoid mojo by keeping the use of WiggleAll to a minimum.

The script performs well in the early puzzle stages, making the backbone more natural and stable. It helps lock in the protein structure right after manual folding, before switching to standard DRW recipes.
Its approach is to generate a large number of solutions (hundreds or even thousands) and evaluate them using the original scoring system without performing a Fuse.

There are two key settings:
1) “Solutions to search”
The number of Rebuild solutions generated in a single run (one pass of the protein).
The default value is twice the number of protein segments, i.e., 100–300 rebuilds.
2) “Make Fuse” (slider):
Controls how often fusing happens. Specifically, it triggers a Fuse after every X iterations that don’t score points.
This is the easiest way to adjust the depth and speed of the Rebuild.

Have fun!
]]--

version = "1.0"

function ScoreReturn()
	x = current.GetEnergyScore()
	return x-x%0.001
end
  -----------------------------------------------------------------------------------------------init
proteinLength=structure.GetCount ()

-- service slots are 5:
-- slot1 is used as basic one
-- slot100 - as the best unfuzed solution
-- slot98 used to store fuzed version of the highest score solution
-- slot99 used as service slot to manage fuze
-- slot97 used as service slot to store lighly fuzed version of the currently fuzed solution
totalSlotsAvailable = 95
maxRebuildCount= proteinLength * 2
slotsToFuze=12
slotsNoGain = 1
shiftFuze = -1
fuzeAfternoGain = 4	
forceChange = -1
shiftNoGain = -1
stopAfter = -1

multipleOverlap = 2

currentFuze = 1

StartRebuild = 13
EndRebuild = 10

reportLevel=3
sphereRadius=8
shakeOnRank=false
convertLoop=true
worstFirst=true
unfuzeFirst=false
fullProtein=true
resetOnGain = false
fuze2=false
--saveLocal=true
remixNotRebuild = false
rarityManagement = false

doFuze = false

endCI=behavior.GetClashImportance()
CI = endCI

energy2BBScoreRatio=2
selectionLength=13

fuzeConfig = {
						"{0.05, -7, 0.05,2} {0.25, 3, 1.00,20} {0.05, 3, 0.25,7} {0.25, 3, 0.25,20} {1.00, 3, 1.00,20}", --best score
						--"{0.25, 1, 1.00,20} {0.05, 2, 0.05,7} {0.25, 3, 0.05,7} {1.00, -2, 1.00,20} {1.00, 3, 1.00,20}", --best score
						"{0.05, -20, 0.05,7} {0.25, 2, 0.25,7} {1.00, 2, 1.00,20}",	--fastest
						"{0.25, -7, 0.05,7} {0.05, 3, 0.25,20} {1.00, 2, 1.00,20}",
						"{1.00, 3, 1.00,7} {0.05, 1, 0.25,2} {1.00, 3, 1.00,20}",
						"{0.05, -2, 0.05,7} {1.00, 3, 1.00,7} {0.05, 2, 0.05,7} {0.25, 1, 0.25,2} {1.00, 3, 1.00,20}"
						}	--fastest

currentFuze = 3

startingAA=1
bestScore=ScoreReturn()
lastScore=-999999
bestFoundScore=-999999
currentBBScore=-999999
medianScore = -999999
medianScore_prev = -999999
rebuildScores = {}
remixBBScores = {}, {}

rebuildScores[1]=bestScore
bestSelectNum=0
remixNum = #rebuildScores
undo.SetUndo(false)


resetProtein=true 
splitDirection = "forward"
bestSlot=0

selectionEnd=0
selectionStart=0
selNum=0
roundsWithNoGain = 0
RoundNumber = 1
totalRebuildsDone = 0
currentSlot = 0
activityMap = {}


selDialog=selNum
selectionStartArr={}
selectionEndArr={}
solutionIDs = {}
action="Rebuild"

numPieces = math.floor(proteinLength/13+0.5)
rebuiltResidues = {} 		--array to monitor changes in residues in all the solutions in slots

startScore=ScoreReturn()
startRoundScore = ScoreReturn()
initScore=ScoreReturn()  --the real start score that doesn't changes

print ("+++Starting score "..startScore.." saved to slot 1")

save.Quicksave(1)
save.Quicksave(98)
save.Quicksave(100)

recentbest.Save()

selection.DeselectAll() --Clear All The Selections On The Start for DRW!
save.SaveSecondaryStructure()
-----------------------------------------------------------------------------------------------main-----------------------------------------------------------------------------------------------
function main()
	resetRebuiltResidues()
	
    FindSelection()
    tempSelNum=selNum
    
    if selNum==0 then selNum=1 end
    selDialog=selNum
    --Call Dialog with the options
    while requestResult~=1 do
		requestResult = RequestOptions()

      if requestResult == 0 then return end -- cancel button pressed by user
      if requestResult == 2 then selDialog=selDialog+1 end
      if requestResult == 3 then 
        SelectLoops() 
        FindSelection()
        --print (selNum,selDialog)
        selDialog=selNum
      end
      if selDialog>5 then selDialog=5 end
    end
    selNum=math.max(tempSelNum, selDialog, selNum)

	--Fuze first if unfuzed solution checkbox is selected
	if (useSlot98) then 
		if reportLevel>2 then print ("Making Fuze for this unfuzed start", startScore) end
		FuzeLocal(fuzeConfig[currentFuze])
		if reportLevel>1 then print ("Fuzed to", ScoreReturn()) end
		if ScoreReturn() > startScore then
			save.Quicksave(98)
			startScore=ScoreReturn()
			startRoundScore = ScoreReturn()
			initScore=ScoreReturn() 
		end
	else
		save.Quicksave(98)
	end
	save.Quickload(1)

    --Report settings in output window
    strOutput=action..": "..maxRebuildCount.." / fuze:"..slotsToFuze..". Range "..StartRebuild.." /"..EndRebuild
    if worstFirst then strOutput=strOutput..", worst first" end
    if resetOnGain then strOutput=strOutput..", reset" end
    if shakeOnRank then strOutput=strOutput..", shake" end
    if reportLevel>2 then strOutput=strOutput..". ReportLevel="..reportLevel end
	strOutput = strOutput.."\n"
    --if useSlot98 then strOutput=strOutput..", useSlot98" end
    if fuze2 then strOutput=strOutput..", fuze2" end
    if slotsNoGain then strOutput=strOutput..", slotsNoGain="..slotsNoGain end
    if forceChange then strOutput=strOutput..", forceChange="..forceChange end
    if fuzeAfternoGain then strOutput=strOutput..", fuzeAfternoGain="..fuzeAfternoGain end
    if shiftFuze then strOutput=strOutput..", shiftFuze="..shiftFuze end
    --if shiftNoGain then strOutput=strOutput..", shiftNoGain="..shiftNoGain end
    if stopAfter>0 then strOutput=strOutput..", stopAfter="..stopAfter end
    if (reportLevel>1) then print (strOutput) end


	print ("-------------------------------------------------")

    loopIter=1
    while (roundsWithNoGain <= stopAfter) or (stopAfter<0) do --for Infinite execution option
		local startTime = os.clock()

		loopIter=loopIter + 1

		if (reportLevel>1) then print("Iteration", (loopIter-1)..".") end

		if (startScore > startRoundScore) then 
			roundsWithNoGain = 0
			ratio_rarity_changed = ratio_rarity --return rarity ratio to initial values
		else
			if loopIter > 2 then --don't change the roundsWithNoGain when just entered the loop
				roundsWithNoGain = roundsWithNoGain + 1
				--increase the part of the rare solutions
				ratio_rarity_changed = math.min (2, ratio_rarity_changed  + 0.2)
				if reportLevel>1 then print ("Score "..startScore," roundsWithNoGain "..roundsWithNoGain) end
			end
		end

		--Make a change in settings every roundsWithNoGain rounds that gained zero points
		if (loopIter>2) then
			-- save Fuzed version of the protein every roundsWithNoGain rounds that gained zero points
			if (roundsWithNoGain >= fuzeAfternoGain)  and (fuzeAfternoGain>0) then
				if reportLevel>1 then print ("NoGain "..roundsWithNoGain.." rounds. fuzeAfternoGain "..fuzeAfternoGain) end
				if startScore ~= initScore then 
					-- if we are stuck then we accept bestslot solution as new and reset all the slots
					lastScore = ScoreReturn()
					save.Quickload(98)
					if ScoreReturn() > lastScore then 
						save.SaveSolution("defuse slot98 "..ScoreBBReturn())
						save.Quicksave(100)
						RoundNumber = RoundNumber + 1
						resetProtein = true
						if reportLevel>1 then print ("Switched to Fuzed version bb ", ScoreBBReturn(), "score ", ScoreReturn()) end
					else 
						save.Quickload(100)
					end
				end
			end

			-- shift between two types of fuze
			if (roundsWithNoGain >= shiftFuze) and (shiftFuze>0) then
				if currentFuze == 2 then 
					currentFuze = 1
					if reportLevel>1 then print ("Switched to Fuze1. (rounds with no gain="..shiftFuze..")") end
				else
					if reportLevel>1 then print ("Switched to Fuze2. (rounds with no gain="..shiftFuze..")") end
					currentFuze = 2
				end
			end
			-- reset every roundsWithNoGain rounds that gained zero points
			if (roundsWithNoGain >= forceChange)  and (forceChange>-1) then
				lastScore = ScoreReturn()
				save.Quickload(100)
				--save.SaveSolution("drw slotXX "..ScoreBBReturn())
				resetProtein = true --reset the slots, so (bestSlot) will be rebuild now.
				if reportLevel > 1 then print("Accepting best from previous Round. NoGain " .. roundsWithNoGain .. " rounds. forceChange " .. forceChange) end
				if reportLevel > 1 then print("Switched to version bb ", ScoreBBReturn(), " score ", ScoreReturn()) end
			end

			if (roundsWithNoGain >= slotsNoGain) and (slotsNoGain>=0)  then
				maxRebuildCount = math.floor(maxRebuildCount * 1.1)
				slotsToFuze = math.floor(slotsToFuze * 1.1) 
				if slotsToFuze > 36 then slotsToFuze = 36 end
				if reportLevel>1 then print ("slotsNoGain reached. Increased maxRebuildCount="..maxRebuildCount..", slotsToFuze="..slotsToFuze) end
			end

			if (roundsWithNoGain >= shiftNoGain) then
				selectionLength = selectionLength - 1  
				if selectionLength < EndRebuild then selectionLength = EndRebuild end 
				if reportLevel>1 then print ("NoGain "..roundsWithNoGain.." rounds. shiftNoGain "..shiftNoGain..". Length is "..selectionLength) end
				--selectionStartArr, selectionEndArr = SplitProteinOnInit(activityMap, numPieces)
				--roundsWithNoGain = 0
				startingAA = 1
			else  
				--shift selections start/end points to prevent rebuilding the same selections over and over 
				--current pointer for the selections markup on full rebuild (shift by one every iteration to prevent rebuilding the same pieces again and again)
				startingAA=startingAA-1 
			end
		end
	
		--------'REBUILD PROTEIN BY SELECTIONS--------------------------------------------------------------------------
		if resetProtein == true then
			for k=1, proteinLength do
				activityMap[k] = 1
			end
		else
			activityMap=getActivityMap()
			printActivity (activityMap)
		end

		--Split protein on selections in 2 diferent ways depending if that is the first iteration
		if loopIter>0 then --first iteration only
			selectionStartArr, selectionEndArr = SplitProteinOnInit(numPieces) --create the selections
		else
			selectionStartArr, selectionEndArr = SplitProteinByActivity(activityMap, numPieces)
		end

		--Calculate number of rebuilds for each selection and adjust maxRebuildCount fo split between sels as int
		--total number of rebuilds remains constant per Iteration, but they are distributed between selection
		tryFirstX, numOfRebuilds, adjustedTotalRebuilds = adjustRebuilParameters(maxRebuildCount, selNum)
		if reportLevel>1 then 
			print ("Rebuilds adjusted to ~"..adjustedTotalRebuilds.." or ~" ..numOfRebuilds.." per sel. " ..tryFirstX.." best slots will be tested.") 
		end
		rebuildsDistribution = DistributeRebuildsByActivity(activityMap, selNum * numOfRebuilds)

		sortedStartArr, sortedEndArr = SortSelections(selectionStartArr, selectionEndArr) --sort array ascending to rebuild worst segments first
		startRoundScore = startScore

		if (reportLevel>2) then printSelections() end
		-------------------
		for m=1, selNum do
			selectionStart=sortedStartArr[m]  --selectionStartArr[m]
			selectionEnd=sortedEndArr[m] 	--selectionEndArr[m]
			
			SetSelection()
			if convertLoop then 
			  save.LoadSecondaryStructure()
			  structure.SetSecondaryStructureSelected("l")
			end
			if resetProtein then 	--make this on first iteration at first selection only, or on solution reset
				resetProtein = false
				resetRebuiltResidues(switchedTo) --reset the monitoring array too(except switchedTo slot)
				resetInheritanceTable(bestSlot)
				solutionsFound=1
				remixBBScores = {}
				rebuildScores = {}
				topXinArray = {}
				currentSlot = 1
			end

			if reportLevel>1 then 
				print((loopIter-1).."."..m.."/"..selNum.." len "..(selectionEnd - selectionStart  +1).." segs: "..selectionStart.."-"..selectionEnd..
						". R"..RoundNumber) --..". Current ScoreBB "..ScoreBBReturn()
			end

			if (#remixBBScores<tryFirstX) then --if not enough solutions yet - just rebuild, dont make that complex as in section lower
				localRebuildsCounter, successRebuilds = RebuildToSlots(rebuildsDistribution[m])
				if reportLevel>1 and reportLevel<4 and localRebuildsCounter>0 then print ("Accepted "..successRebuilds.." solutions")  end
				if reportLevel > 3 then print ("--------------------")	end
			else
				topXinArray = returnXinOrder(remixBBScores, tryFirstX)
				save.Quickload(topXinArray[1])
				activityMap=getActivityMap()
				if reportLevel > 3 then printActivity(activityMap) end

				--now for each of the best X slots try creating best rebuild solutions for current selection
				for currSel=1, tryFirstX do
					--in case when tryFirstX > existed Rebuild slots, we prevent the recipe crash on the first iteration pass
					if currSel > #topXinArray then 
						if reportLevel>2 then print (#topXinArray.." slots are not enoug. Adding some more ++++") end
						topXinArray = returnXinOrder(remixBBScores, tryFirstX) end 
					currentSlot = topXinArray[currSel] --the number of slot that is now rebuilding

					save.Quickload(currentSlot)
					successRebuilds = RebuildToSlots(rebuildsDistribution[m])
					if reportLevel>2 then print ((loopIter-1).."."..m.."/"..selNum.." "..currSel.."/"..tryFirstX.." ",
						"Slot "..currentSlot, " accepted ", successRebuilds, " / "..rebuildsDistribution[m], "bb"..ScoreBBReturn()) end
					if reportLevel > 3 then print ("----------------------------------------") end
				end
			end
			activityMap=getActivityMap()

			if reportLevel > 3 then printActivity (activityMap) end
		end --for selNum


		bestScore=-999999
		bestSlot=0
		bestSelectNum = math.min (slotsToFuze, #remixBBScores)
		activityMap=getActivityMap()
		if reportLevel > 3 then printActivity(activityMap) end
		if reportLevel>1 then print("Total Rebuilds Tested:", totalRebuildsDone, "Ranking best "..bestSelectNum) end

		SortByBackbone()
		printHistory()
		local endTime = os.clock()
		executionTime = endTime - startTime
		if reportLevel>2 then print (string.format("----------------------%.2f %.2f sec.-------------------", executionTime, executionTime2)) end
		executionTime = 0
		executionTime3 = 0
		executionTime3 = 0
		
		--------'FUZE BEST SOLUTIONS
		if #remixBBScores > 0 then 

			lastSlots = {} --to restore best slots on clean up
			shortFuzeScores = {} -- to keep track of scores after short fuze

			--Short Fuze a lot of slots and score best
			testedSlots = 0
			for i, remixBBScore in ipairs(remixBBScores) do
				testedSlots = testedSlots + 1
				if testedSlots > slotsToFuze then break end

				--Test only non-fuzed slots
				if  remixBBScore.score ~= 0 then  
					testedSlots = testedSlots - 1 
				else			
					table.insert(lastSlots, remixBBScore.id)
					save.Quickload(remixBBScore.id)
					prefuzeBB = ScoreBBReturn()

					FuzeLocal(fuzeConfig[2]) -- short fuze

					currentScore = ScoreReturn()
					currentBBScore = ScoreBBReturn()
					remixBBScore.score = currentScore --update remixBB array values with the fuzed ones
					descriminator = currentBBScore - prefuzeBB
					UpdateRemixBBDescr(i, descriminator)
					updateHistoryArray (remixBBScore.id, currentScore, currentBBScore)
					
					improve = roundX(currentScore - startScore)
					textBest = ""
					if currentScore > bestScore  then 
						save.Quicksave(97) 
						save.Quickload(remixBBScore.id) --delete?
						bestScore = currentScore
						bestSlot = remixBBScore.id
						textBest = "*"
					end

					table.insert(shortFuzeScores, {id = remixBBScore.id, score = currentScore})

					if reportLevel > 1 then 
						print("Fuzed slot "..remixBBScore.id, " from bb "..prefuzeBB.." to bb "..currentBBScore, "score: "..currentScore, "/"..remixBBScore.rank, textBest)
					end
					if reportLevel > 3 then printRow(remixBBScore.id) end
				end
			end

			-- This section is to fuze best results with some better fuze
			if  (#shortFuzeScores>0) then 
				-- Long Fuze sqrt of best short fuzes and find the best one
				table.sort(shortFuzeScores, function(a, b) return a.score > b.score end)
				local numToLongFuze = math.floor(math.max(1, math.sqrt(slotsToFuze)))  --switch min to max
				numToLongFuze = math.min(numToLongFuze, #shortFuzeScores) --prevent crash when no unfuzed left

				for i = 1, numToLongFuze do
					local bestShortFuze = shortFuzeScores[i]
					save.Quickload(bestShortFuze.id)

					FuzeLocal(fuzeConfig[1]) -- long fuze

					currentScore = ScoreReturn()
					currentBBScore = ScoreBBReturn()
					if currentScore > bestShortFuze.score and false then --disabled
						--update remixBB and History array values
						local idxx = getRemixBBIDBySlot(bestShortFuze.id)
						remixBBScores[idxx].score = currentScore 
						descriminator = currentBBScore - prefuzeBB --descriminator<0 means we have some disassembly
						UpdateRemixBBDescr(idxx, descriminator)
						updateHistoryArray (remixBBScore.id, currentScore, currentBBScore)
					end

					if currentScore > bestScore then 
						bestScore = currentScore
						bestSlot = bestShortFuze.id
						save.Quicksave(97)
						if reportLevel > 1 then 
							print("Best slot after long fuze "..bestShortFuze.id, " score: "..currentScore)
						end
					end
				end --for

				------------ Manage best slot
				save.Quickload(97)
				if reportLevel>1 then print ("Fuzed slot "..bestSlot.." bb "..ScoreBBReturn(), "score ",ScoreReturn() ) end

				if reportLevel>2 then printScoreStats() end	
				--print ("Residues that are not rebuilt are shown with 0:")
				temp = getRow(bestSlot)
				if reportLevel>2 then printRow(bestSlot) end --report what residues had been changed in this slot

				lastScore = ScoreReturn()
				lastScoreBB = ScoreBBReturn()

				if lastScore > bestFoundScore then 
					--save.SaveSolution("defuze best slot"..bestSlot.."  "..lastScore)
					if reportLevel>2 then print ("saved best found slot"..bestSlot.." score "..lastScore) end
				end

				if (ScoreReturn() > startScore) then
					save.Quicksave(98) 
					if reportLevel>1 then print ("Gained "..roundX(ScoreReturn()-startScore).." points. New score is", ScoreReturn()) end
					--if reportLevel>1 then print ("Replaced old "..startScore..". Score improved to bb "..ScoreBBReturn(), "score: ", ScoreReturn()) end
					--save.SaveSolution("defuse slot"..bestSlot.." bb "..lastScoreBB)
					if reportLevel>3 then print ("saved drw slot"..bestSlot.." bb "..lastScoreBB) end
					startScore = ScoreReturn()
					--resetProtein = true
					switchedTo = bestSlot
					save.Quickload(bestSlot)
				else
					save.Quickload(100)
					--if reportLevel>2 then print ("Best score "..startScore..". New one is not enough "..ScoreReturn()) end
					if reportLevel>2 then print ("No gain, Restoring "..startScore ) end
				end
				
				if lastScore > bestFoundScore then 
					save.Quicksave(100) --save unfuzed version of the best solution found
					bestFoundScore = lastScore
				end

				--print ("Best score "..startScore.." Fuze score /bb "..lastScore.." "..lastScoreBB, "preFuze score /bb "..ScoreReturn(), ScoreBBReturn() )
			end  --if  (#shortFuzeScores>0) 
		end --if #remixBBScores > 0
		if reportLevel>1 then print (string.format("----------------------%.2f sec.-------------------", executionTime3)) end
    end --while

	print ("Initial version and fuzed one are in undo stack")
    
    Cleanup()
end -- function main()

------------------------------------------------------------------------------------------------Dialog-----------------------------------------------------------------------------------------------

function RequestOptions()
    ask=dialog.CreateDialog("Options: Fast Assembly"..version)

    ask.maxRebuildCount = dialog.AddSlider("Solutions to search",maxRebuildCount,1,15000,0) --up to 98 slots possible (change 36 to 98 if needed)
    ask.slotsToFuze = dialog.AddSlider("Slots to fuze",slotsToFuze,1,36,0)

    ask.remixNotRebuild=dialog.AddCheckbox("Remix instead of Rebuild", remixNotRebuild)
    --ask.remixNotRebuild.value=false
    ask.convertLoop=dialog.AddCheckbox("Convert to Loop", convertLoop)
    ask.resetOnGain=dialog.AddCheckbox("Reset slots on Gain", resetOnGain)

    ask.numPieces = dialog.AddSlider("numPieces",numPieces,2,math.floor(proteinLength/3),0)
    ask.StartRebuild = dialog.AddSlider("Start Length",StartRebuild,2,proteinLength,0)
    ask.EndRebuild = dialog.AddSlider("End Length",EndRebuild,2,proteinLength,0)
    ask.multipleOverlap = dialog.AddSlider("Overlap",multipleOverlap,0,10,0)

    ask.l1 = dialog.AddLabel("Fuze:")
    ask.shakeOnRank=dialog.AddCheckbox("Shake solution on Rank stage", shakeOnRank)
    ask.worstFirst=dialog.AddCheckbox("Worst First", worstFirst) 
    ask.doFuze=dialog.AddCheckbox("Store Fuzed", doFuze)
    --ask.forceChange=dialog.AddCheckbox("Force Changes (loss)", forceChange)

	ask.l2 = dialog.AddLabel("Do after X rounds with no Gain")
    ask.slotsNoGain = dialog.AddSlider("Slots increase",slotsNoGain,-1,15,0)
	ask.shiftFuze = dialog.AddSlider("Fuze local",shiftFuze,-1,15,0) 	--: change rebuild length when there are shiftNoGain iterations with no gain
	ask.fuzeAfternoGain = dialog.AddSlider("Make Fuze",fuzeAfternoGain,-1,15,0) 	--: change rebuild length when there are shiftNoGain iterations with no gain
	ask.forceChange = dialog.AddSlider("Force change (loss)",forceChange, -1,20,0)  --
	ask.shiftNoGain = dialog.AddSlider("Rebuild length",shiftNoGain,-1,15,0) 	--: change rebuild length when there are shiftNoGain iterations with no gain
	ask.stopAfter = dialog.AddSlider("Stop After",stopAfter,-1,100,0) 	--: change rebuild length when there are shiftNoGain iterations with no gain

    ask.fuze2=dialog.AddCheckbox("Fuze2", false)
    ask.useSlot98=dialog.AddCheckbox("Unfuzed version", false)
    ask.reportLevel = dialog.AddSlider("Report detalization", reportLevel,1,5,0)
 
    ask.OK = dialog.AddButton("OK",1) 
    --ask.addSelections = dialog.AddButton("AddSelection",2) 
    ask.selectLoops = dialog.AddButton("SelLoops",3) 
    ask.Cancel = dialog.AddButton("Cancel",0)

    returnVal=dialog.Show(ask)
	if returnVal > 0 then
		
		if returnVal==1 then maxRebuildCount=ask.maxRebuildCount.value end

		remixNotRebuild=ask.remixNotRebuild.value
		worstFirst=ask.worstFirst.value
		shakeOnRank=ask.shakeOnRank.value
		reportLevel=ask.reportLevel.value
		resetOnGain=ask.resetOnGain.value
		convertLoop = ask.convertLoop.value
		slotsToFuze=ask.slotsToFuze.value
		forceChange=ask.forceChange.value
		doFuze=ask.doFuze.value
		fuze2=ask.fuze2.value
		stopAfter=ask.stopAfter.value
		slotsNoGain=ask.slotsNoGain.value
		fuzeAfternoGain=ask.fuzeAfternoGain.value
		useSlot98 =ask.useSlot98.value

		multipleOverlap = ask.multipleOverlap.value
		forceChange=ask.forceChange.value 
		shiftNoGain=ask.shiftNoGain.value 
		shiftFuze=ask.shiftFuze.value 

		numPieces=ask.numPieces.value
		StartRebuild=ask.StartRebuild.value
		EndRebuild=ask.EndRebuild.value

		----fix some vars if needed
		if remixNotRebuild then 
			action="Remix" 
			if StartRebuild > 9  then 
				StartRebuild = 9 
				print ("StartRebuild more than 8 not available for Remix. Setting at 8")
			end
			if EndRebuild > 9  then 
				EndRebuild = 9 
				print ("EndRebuild more than 8 not available for Remix. Setting at 8") 
			end
		end

		if StartRebuild < EndRebuild then 
			temp = StartRebuild
			StartRebuild = EndRebuild
			EndRebuild = temp
		end
		selectionLength=StartRebuild

		if slotsToFuze > maxRebuildCount then slotsToFuze = maxRebuildCount end

		selectionStartArr, selectionEndArr = SplitProteinOnInit(numPieces) --create the selections

		-- to add 0.5 to every highest subscore solution (round first decimal "bestSelectNum / 12". 10 is used as the first decimal floor basis)

    else 
      print ("Canceled") 
    end

	return returnVal
end


---------------------------------------------------------------------------------------------Rebuild/Remix---------------------------------------------------------------------------------------------
generation = 0
function RebuildToSlots(Slots2Rebuild)
	generation = generation + 1
	remixNum=#rebuildScores-1
	j=0
	rebuildIter=1
	repeats_counter=0
	undo.SetUndo(false) 
	localRebuildsCounter = 0

	save.Quicksave(99)
	if reportLevel>3 then print ("Rebuilding slots:",Slots2Rebuild) end
	while (j<Slots2Rebuild) and (rebuildIter<7) do

		j=j+1
		structure.RebuildSelected(rebuildIter)

		isDuplicated = CheckRepeats() 
		if isDuplicated then 
			j=j-1
			rebuildIter = rebuildIter+1
			repeats_counter=repeats_counter+1

		else
			totalRebuildsDone = totalRebuildsDone + 1
			lastScore = ScoreReturn()
			if shakeOnRank then TinyFuze() end --makes a small shake before ranking energy score of rebuild (if checkbox was selected).
			currentRow = getRow(currentSlot)
			currentRow = changeRowValues(currentRow, selectionStart, selectionEnd)
			currentScore=ScoreReturn()
			currentScoreBB=ScoreBBReturn()
			
			parentID =  getRemixBBIDBySlot(currentSlot)
			descriminator = 0
			if parentID~=0 and parentID~=nil then 
				descriminator = remixBBScores[parentID].bbDescriminator 
			end

			--Two different scenarios: 
			--1) when the solution array is empty (first rebuilds just adds values to arrays), 
			--2) and when it is full (then replace the min value instead of inserting)
			if #rebuildScores <= totalSlotsAvailable then
				localRebuildsCounter = localRebuildsCounter + 1
				solutionsFound = solutionsFound+1
				rebuildScores[solutionsFound] = lastScore --saving score before TinyFuze to really check the repeats for next rebuilds
				save.Quicksave(solutionsFound)
				subScore = GetSolutionSubscores() --saving all the subscores of current solution in an Array to find later the one with the highest value on each subscore
				local newEntry = {id = solutionsFound, scoreBB = currentScoreBB, bbDescriminator=descriminator, score = 0, subScores = subScore,  generation=generation,  rank = 0, rarity = 0}
				table.insert(remixBBScores, newEntry)
				
				--printActivity (currentRow)
				addRow(currentRow)
				addHistoryEntry(solutionsFound, currentSlot, currentScoreBB, currentRow, generation)

				if reportLevel>3 then print (j..". Slot", solutionsFound, "backbone", currentScoreBB,  "score", currentScore) end

			else
				--Add every solution to remixBBScores. Then sort and remove the worse one. 
				newElement= {id = 0, scoreBB=currentScoreBB, bbDescriminator=descriminator,score=0, subScores = GetSolutionSubscores(), generation=generation, rank = 0, rarity = 0} 
				table.insert(remixBBScores, newElement)
				SortByBackbone(rarityManagement) --sort by rank
				
				IdxToRemove = #remixBBScores
				minSlotNo = remixBBScores[#remixBBScores].id
				-- Do not add too many from one generation, - replace the worst record only among the same generation
				if  j>tryFirstX then
					for i = #remixBBScores, 1, -1 do
						if remixBBScores[i].generation == generation then
							IdxToRemove = i
							minSlotNo = remixBBScores[i].id
							--localRebuildsCounter = localRebuildsCounter - 1
							break
						end
					end
				end
				
				--if new solution is good enough then get id of the worst
				if remixBBScores[IdxToRemove].id ~= 0 then 
					newElement.id = minSlotNo --it is a reference to remixBBScores last added
					if reportLevel>4 then print ("minSlot, oldBB, oldRank, currBB", minSlotNo, IdxToRemove, remixBBScores[#remixBBScores].scoreBB, remixBBScores[#remixBBScores].rank, currentScoreBB) end

					localRebuildsCounter = localRebuildsCounter + 1
					save.Quicksave(minSlotNo)
					rebuildScores[minSlotNo] = lastScore --saving score before TinyFuze to really check the repeats for next rebuilds
					removeSlotStats(minSlotNo) --remove worst from history table
					addHistoryEntry(minSlotNo, currentSlot, currentScoreBB, currentRow, generation)  --add temporary to calc rarity
					setRow(currentRow, minSlotNo)
					if reportLevel>4 then printRow(minSlotNo) end
					if reportLevel>3 then print (j..". Replaced Slot "..minSlotNo.." bb "..currentScoreBB.." score "..currentScore) end
				end
				table.remove(remixBBScores, IdxToRemove)
			end
			rebuildIter=1
		end
		--clean history table periodicaly to mantain the calculation speed
		if totalRebuildsDone % 300 ==  0 then cleanHistory() end
		save.Quickload(99)
	end --while
	if reportLevel>3 and (repeats_counter > 5) then print ("repeats_counter",repeats_counter) end
	undo.SetUndo(true) 
return localRebuildsCounter, j
end

executionTime2 = 0
executionTime3 = 0

ratio_scoreBB = 0.5
ratio_score = 1
ratio_subScores = 0.5
ratio_rarity = 0.2
ratio_rarity_changed = ratio_rarity

function SortByBackbone(calcRarity)
    local startTime2 = os.clock()

	local bestSelectNum = #remixBBScores --rank points value 


	for i, remixBBScores in ipairs(remixBBScores) do
		remixBBScores.rank=0  --clear array
	end
	--rank best backbone score
	table.sort(remixBBScores, function(a,b) return (a.scoreBB + a.bbDescriminator)< (b.scoreBB + b.bbDescriminator) end)
	for i, remixBBScores in ipairs(remixBBScores) do
		remixBBScores.rank=remixBBScores.rank + i*ratio_scoreBB 
	end
	--rank best energy score 
	-- Secondary sort by scoreBB if scores are equal (score==0)
	table.sort(remixBBScores, function(a, b)
		if a.score == b.score then	return a.scoreBB < b.scoreBB 
		else return a.score < b.score end
	end)
	for i, remixBBScores in ipairs(remixBBScores) do
		remixBBScores.rank=remixBBScores.rank + i*ratio_score
	end	
	--rank rarity
	if calcRarity~=nil  and calcRarity==true then
		--fill remixBBScores with rarity values
		for i, remixBBScores in ipairs(remixBBScores) do
			local historyIdx = historyIDbySlot[remixBBScores.id]
			--print (remixBBScores.id, historyIdx)
			remixBBScores.rarity = calculateRarityScore(historyIdx)
		end
		--rank by rarity
		table.sort(remixBBScores, function(a,b) return a.rarity < b.rarity end)
		for i, remixBBScores in ipairs(remixBBScores) do
			remixBBScores.rank=remixBBScores.rank + i*ratio_rarity_changed
		end
	end
	-- rank for solutions that have one of top subscores
	highestSolutionIDs, subScoresNumber = GetHighestSolutionIDs()
    for _, remixBBScores in ipairs(remixBBScores) do
        local solutionID = remixBBScores.id
        for _, solutionIDs in pairs(highestSolutionIDs) do
            for _, id in ipairs(solutionIDs) do
                if id == solutionID then
                    remixBBScores.rank = remixBBScores.rank + bestSelectNum/subScoresNumber * ratio_subScores	
                end
            end
        end
    end
	table.sort(remixBBScores, function(a,b) return a.rank > b.rank end)
	for i, remixBBScores in ipairs(remixBBScores) do
		remixBBScores.rank=math.floor(remixBBScores.rank * 10 + 0.5) / 10 
	end

	if (reportLevel>4) and (calcRarity~=nil) then 
		for scorePart, solutionIDs in pairs(highestSolutionIDs) do
			print("Solution ID with highest",scorePart," subscore:", table.concat(solutionIDs, ", "), subScoresNumber)
		end
		for i, remixBBScores in ipairs(remixBBScores) do
			print("Slot",remixBBScores.id,"bb:", remixBBScores.scoreBB, "score", remixBBScores.score)
		end
	end
    local endTime2 = os.clock()
    executionTime2 = executionTime2 + endTime2 - startTime2
end

-- Function to get the list of SolutionIDs with the highest subscores for each score part
function GetHighestSolutionIDs()
    local highestSolutionIDs = {}
    local solutionSubscoresArray = {}
	local counterx = 0
    for _, entry in ipairs(remixBBScores) do
        table.insert(solutionSubscoresArray, {id = entry.id, subScores = entry.subScores})
    end

    -- Iterate over each score parts
    for scorePart, _ in pairs(solutionSubscoresArray[1].subScores) do
		local maxSubscore = -999999
		local minSubscore =  999999
		local maxSolutionIDs = {}

		-- Iterate over each solution subscores
		for _, subscores in ipairs(solutionSubscoresArray) do
			local subscore = subscores.subScores[scorePart]
			local solutionID = subscores.id
			-- find the max value between all the solutions for secific subscore
			if subscore > maxSubscore then
				maxSubscore = subscore
				maxSolutionIDs = {solutionID}
			elseif subscore == maxSubscore then
				table.insert(maxSolutionIDs, solutionID)
			end

			-- Track the minimum subscore to prevent including in highestSolutionIDs array subscores with the same values for all the solutions
			if subscore < minSubscore then
				minSubscore = subscore
			end
        end

        if maxSubscore ~= minSubscore then
			highestSolutionIDs[scorePart] = maxSolutionIDs
			counterx = counterx + #maxSolutionIDs
		end
    end

    -- Return the list of SolutionIDs with the highest non-zero subscore for each score part
    return highestSolutionIDs, counterx
end


function CheckRepeats()
	currentScore=ScoreReturn()
	isDuplicate=false
	remixNum=#rebuildScores
	  
	for k=1, remixNum do
		if rebuildScores[k] == currentScore then 
			isDuplicate=true 
			if reportLevel>3 then print (currentScore, "is duplicated to already found in slot",k) end
		end
	end
	return isDuplicate
end

--set Descriminator for current element and all its children in remixBBScores array
function UpdateRemixBBDescr(idxRemixBB, descriminator,level)
	--print (idxRemixBB, descriminator)
	slot = remixBBScores[idxRemixBB].id
	if level==nil then level=0 end
	if reportLevel>4 then
		local indent = string.rep("\t", level)
		print(indent..slot, descriminator)
	end
	remixBBScores[idxRemixBB].bbDescriminator = descriminator
	CurrentID = historyIDbySlot[slot]
	 --change descriminator for non-scored records only
	if (history[CurrentID].score == 0) then 
		for _, offspringID in ipairs(history[CurrentID].offspring) do
			if offspringID and history[offspringID] and history[offspringID].slot then
				slot = history[offspringID].slot
				if (slot ~= 0) then --if the record is active
					idxRemixBB = getRemixBBIDBySlot(slot)
					if idxRemixBB==nil then print ("-", offspringID, slot, CurrentID) end
					UpdateRemixBBDescr(idxRemixBB, descriminator, level+1)
				end
			else
				print("Error: history, offspringID, or slot is invalid")
			end
		end
	end
end

-- Function to get an entry by its ID
function getRemixBBIDBySlot(slot)
	if slot == 1 then return 0 end
	for idx, remixBBScore in ipairs(remixBBScores) do
		if remixBBScore.id == slot then
			return idx
		end
	end
    return nil -- Return nil if no entry with the given ID is found
end
function getHistoryIDBySlot(slot)
	for idx, info in ipairs(history) do
		if info.slot == slot then
			return  info.id
		end
	end
    return nil -- Return nil if no entry with the given ID is found
end


----------------------------------------------Rebuild parameters
function adjustRebuilParameters(totalRebuildCount, selNum)
    local numOfRebuilds
    local tryFirstX

    -- Calculate numOfRebuilds
    if totalRebuildCount / selNum <= 27000 then
        -- Case where sqrt(numOfRebuilds) <= 30
        numOfRebuilds = (totalRebuildCount / selNum) ^ (2/3)
        tryFirstX = math.min(math.sqrt(numOfRebuilds), 30)
    else
        -- Case where sqrt(numOfRebuilds) > 30
        numOfRebuilds = math.floor(totalRebuildCount / (30 * selNum))
        tryFirstX = 30
    end
	tryFirstX, numOfRebuilds = math.floor(tryFirstX + 0.5), math.floor(numOfRebuilds + 0.5)
    local recalculatedTotal = tryFirstX * numOfRebuilds * selNum

    return tryFirstX, numOfRebuilds, recalculatedTotal
end
--------------------------------------------------Set of service functions to monitor what parts of the protein in the accepted solution where actually rebuilt
--Using 2d array. Number of rows = number of occupied slots
--Number of columns = number residues. Each row contain zero value for the residues that werent changed and 1 if they were rebuild in this solution

-- Function to add a new row with'1' between startIndex and endIndex and 0 elsewhere
function addRow(newRow)
    local numColumns = rebuiltResidues.numColumns or 0
    -- Create a row with zeroes if newRow is not provided
    if newRow == nil then
        newRow = {}
        for i = 1, numColumns do
            table.insert(newRow, 0)
        end
    end
    table.insert(rebuiltResidues, newRow)
end
-- Function to modify the row and increase values between selectionStart and selectionEnd by 1
function changeRowValues(row, selectionStart, selectionEnd)
	for j = selectionStart, selectionEnd do
		row[j] = row[j] + 1
    end
    return row
end

-- Function to modify the row and increase values between selectionStart and selectionEnd by 1
-- If 9 is reached, continue with 'a', 'b', ..., 'z', then 'A', 'B', ..., 'Z'
function changeRowValues2(row, selectionStart, selectionEnd)
    for j = selectionStart, selectionEnd do
        if j >= 1 and j <= #row then
            if row[j] == 9 then
                row[j] = 'a'
            elseif type(row[j]) == 'string' then
                if row[j] >= 'a' and row[j] < 'z' then
                    row[j] = string.char(string.byte(row[j]) + 1)
                elseif row[j] == 'z' then
                    row[j] = 'A'
                elseif row[j] >= 'A' and row[j] < 'Z' then
                    row[j] = string.char(string.byte(row[j]) + 1)
                elseif row[j] == 'Z' then
                    row[j] = 'Z'
                end
            elseif type(row[j]) == 'number' and row[j] >= 0 and row[j] < 9 then
                row[j] = row[j] + 1
            end
        end
    end
    return row
end

-- Function to perform a shallow copy of a table (array) and return the_copy and not the_link to the original row from rebuiltResidues[]  in getRow()
local function shallowCopy(array)
    local copy = {}
    for i, v in ipairs(array) do
        copy[i] = v
    end
    return copy
end
-- Function to get a copy of a row
function getRow(rowNumber)
    if rebuiltResidues[rowNumber] then
        return shallowCopy(rebuiltResidues[rowNumber])
    end
end
function setRow(sourceRow, destRowNumber)
    if rebuiltResidues[destRowNumber] then
        for j = 1, #sourceRow do
            rebuiltResidues[destRowNumber][j] = sourceRow[j]
        end
    end
end
-- Function to print the values of the specified row compactly

function resetRebuiltResidues(rowNumber)
	if (reportLevel>1) then print ("Stats reset") end
	if rowNumber~=nil and rowNumber~=0 then
		currentrow = getRow(rowNumber)	--save the row before reseting the table
	end
	rebuiltResidues = {} 
	rebuiltResidues.numColumns = proteinLength
	addRow() --add and empty row for 'slot1'

	if rowNumber~=nil  and rowNumber~=0 then
		setRow(currentrow, 1)
		print (table.concat(rebuiltResidues[1], ""))
	end
end

function printRow(rowNumber) 
    if rebuiltResidues[rowNumber] then
        print(table.concat(rebuiltResidues[rowNumber], ""))
    end
end
----------------------------------------------------------activityMap stats for the whole rebuild
-- Convertion between chars and number is activity map
function charToValue(c)
    local ascii = string.byte(c)
    if ascii >= string.byte('0') and ascii <= string.byte('9') then
        return ascii - string.byte('0')
    elseif ascii >= string.byte('a') and ascii <= string.byte('z') then
        return ascii - string.byte('a') + 10
    elseif ascii >= string.byte('A') and ascii <= string.byte('Z') then
        return ascii - string.byte('A') + 36
    end
end
function valueToChar(value)
    if value >= 0 and value <= 9 then
        return string.char(string.byte('0') + value)
    elseif value >= 10 and value <= 35 then
        return string.char(string.byte('a') + value - 10)
    elseif value >= 36 and value <= 61 then
        return string.char(string.byte('A') + value - 36)
    else
		if (reportLevel>2) then print ("ActivityMap element out of Range!!", value) end
		return string.char(string.byte('Z'))
	end
end
function normalizeScores(activityScores)
    local minScore = math.huge
    local maxScore = -math.huge

    for _, score in ipairs(activityScores) do
        if score < minScore then minScore = score end
        if score > maxScore then maxScore = score end
    end
	local multiplier = 61 --math.min(61, maxScore) -- 61 iz max letter 'Z'
	--print ("minScore, maxScore",minScore, maxScore, multiplier)

    -- Normalize scores to fit in the range 0 to 61
    local normalizedScores = {}
	for i, score in ipairs(activityScores) do
		normalizedScores[i] = math.floor(((score - minScore) / (maxScore - minScore))*multiplier + 0.5) -- round to the nearest integer
	end
    return normalizedScores
end
function printActivity(activityMap)
	activityMapNormalized = normalizeScores(activityMap)
	activityString = scoresToCharString(activityMapNormalized)
	if (reportLevel>3) and (selNum>0) then print("Activity map:") end 
	if (reportLevel>2) and (selNum>0) then print(activityString) end
end

function scoresToCharString(normalizedScores)
    local result = {}
    for _, score in ipairs(normalizedScores) do
        table.insert(result, valueToChar(score))
    end
    return table.concat(result)
end

-- Function to count top percent results
function getActivityMap()
    local numElements = #rebuiltResidues[1]
    local numTests = #rebuiltResidues-1
    local activityScores = {}
	decayFactor = 0.1
    -- Initialize activityScores
    for i = 1, numElements do
        activityScores[i] = 0
    end
    -- Sum the activity scores weighted by exponential rank
    for testIndex = 1, numTests do
        local rank = remixBBScores[testIndex].rank
        local weight = math.exp(-decayFactor * (rank - 1))
        for elementIndex = 1, numElements do
            local value = rebuiltResidues[testIndex][elementIndex]
            activityScores[elementIndex] = activityScores[elementIndex] + value * weight
        end
    end
    -- Calculate the average weighted score
	local maxValue = 0
    for elementIndex = 1, numElements do
        activityScores[elementIndex] = activityScores[elementIndex] / numTests
		if activityScores[elementIndex] > maxValue then maxValue = activityScores[elementIndex] end
    end
	activityScores = normalizeScores(activityScores)
    return activityScores
end

--------------------------------------------------Other service functions--------------------------------------------------
-- Function to analyze scoreBB in the remixBBScores array
function printScoreStats()
	if #remixBBScores == 0 then return nil end
	
	local highestScore = remixBBScores[1].scoreBB
	local lowestScore = remixBBScores[1].scoreBB
	local sum = 0
	local scores = {}

	for _, entry in ipairs(remixBBScores) do
		local scoreBB = entry.scoreBB
		if scoreBB > highestScore then highestScore = scoreBB end
		if scoreBB < lowestScore then lowestScore = scoreBB end
		sum = sum + scoreBB
		table.insert(scores, scoreBB)
	end

	meanScore = sum / #remixBBScores
	table.sort(scores)
	--local medianScore

	if #scores % 2 == 0 then
		medianScore = (scores[#scores / 2] + scores[#scores / 2 + 1]) / 2
	else
		medianScore = scores[math.ceil(#scores / 2)]
	end

	print("Highest bb "..highestScore.." lowest "..lowestScore..". Mean "..roundX(meanScore, 0).." median "..medianScore)
end

------------------------------------------------------------------------------------------Split protein------------------------------------------------------------------------------------------
--Split protein based on the Activity map. The low acitivty segments are the points where the split is most likely to happen
-- Function to calculate the split points
function calculateSplitPoints(activityMap, numPieces)
    local activityScores = {}

    for i = 1, #activityMap do
        table.insert(activityScores, charToValue(activityMap:sub(i, i)))
    end
    local splitPoints = {}
    local threshold = 20 -- Adjust threshold as needed
    local constant = 30 -- Adjust constant as needed
    for i = 1, #activityScores - 1 do
        local randomFactor = math.random()
        local activity = activityScores[i]
		--activity is about 0..61. so approximately in 64% of cases the split will be accepted.
        if randomFactor * (activity + constant) > threshold then
            table.insert(splitPoints, i)
        end
    end

    -- Ensure that we have exactly numPieces - 1 split points
    while #splitPoints > numPieces - 1 do
        table.remove(splitPoints, math.random(#splitPoints))
    end

    while #splitPoints < numPieces - 1 do
        local candidate = math.random(2, #activityScores - 2)
        if not table.contains(splitPoints, candidate) then
            table.insert(splitPoints, candidate)
        end
    end
    table.sort(splitPoints)
    -- Ensure that each piece has at least 3 elements
    for i = 1, #splitPoints do
        if i == 1 then
            if splitPoints[i] < 3 then
                splitPoints[i] = 3
            end
        else
            if splitPoints[i] - splitPoints[i - 1] < 3 then
                splitPoints[i] = splitPoints[i - 1] + 3
            end
        end
    end
    -- Adjust last split point to ensure the last piece has at least 3 elements
    if #activityScores - splitPoints[#splitPoints] < 3 then
        splitPoints[#splitPoints] = #activityScores - 3
    end
    return splitPoints
end
-- Function to split the structure based on split points
function splitStructure(structure, splitPoints)
    local pieces = {}
    local startIdx = 1

    for _, splitPoint in ipairs(splitPoints) do
        table.insert(pieces, structure:sub(startIdx, splitPoint))
        startIdx = splitPoint + 1
    end

    table.insert(pieces, structure:sub(startIdx))

    return pieces
end
function sortSelections(selectionStartArr, selectionEndArr)
    local combined = {}
    for i = 1, #selectionStartArr do
        table.insert(combined, {start = selectionStartArr[i], finish = selectionEndArr[i]})
    end
    table.sort(combined, function(a, b)
        if a.start == b.start then return a.finish < b.finish
        else return a.start < b.start
        end
    end)
    for i = 1, #combined do
        selectionStartArr[i] = combined[i].start
        selectionEndArr[i] = combined[i].finish
    end
end

-- Main function to split protein by selections
function SplitProteinByActivity(activityMap, numPieces)
	selectionStartArr = {}
	selectionEndArr = {}

	for i = 1, multipleOverlap do
		local splitPoints = calculateSplitPoints(activityMap, numPieces)
		local startIdx = 1
		for _, splitPoint in ipairs(splitPoints) do
			table.insert(selectionStartArr, startIdx)
			table.insert(selectionEndArr, splitPoint)
			startIdx = splitPoint + 1
		end
		table.insert(selectionStartArr, startIdx)
		table.insert(selectionEndArr, #activityMap)
	end

	sortSelections(selectionStartArr, selectionEndArr)
	selNum = #selectionStartArr
    return selectionStartArr, selectionEndArr
end

--------------------Initial protein split on selections
-- Box-Muller transform to generate normally distributed values to form random selLength
function randomGaussian(mean, stddev)
    local u1 = math.random()
    local u2 = math.random()
    local z0 = math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2)
    return z0 * stddev + mean
end
-- Function to generate a random value where the mean is X and 95% of values are within 30% of X
function randomValueWithConfidence(X)
    local mean = X
    local stddev = 0.15 * X
    local value
    repeat
        value = randomGaussian(mean, stddev)
    until value >= 3  -- Ensure the selection length is not less than 3
    return math.floor(value)
end

-- Function to split the protein based on random selection lengths and overlap
function SplitProteinOnInit(numPieces)
    local totalLength, selectionLength = proteinLength, math.floor(proteinLength / numPieces)
    local overlapFraction = 0.4
    local selectionStartArr, selectionEndArr, startIdx, direction = {}, {}, 1, 1

    -- Determine the direction based on the global variable
    if splitDirection == "backward" then
        startIdx, direction = totalLength, -1
        splitDirection = "forward"
    else
        splitDirection = "backward"
    end

    selNum = 0
    while (direction == 1 and startIdx <= totalLength) or (direction == -1 and startIdx >= 1) do
        selNum = selNum + 1
        local currentSelectionLength = randomValueWithConfidence(selectionLength)
        local endIdx = startIdx + (currentSelectionLength - 1) * direction
        endIdx = math.max(1, math.min(totalLength, endIdx))
        if direction == 1 and totalLength - endIdx < 2 then endIdx = totalLength end
        if direction == -1 and endIdx < 3 then endIdx = 1 end
        table.insert(selectionStartArr, (direction == 1 and startIdx or endIdx))
        table.insert(selectionEndArr, (direction == 1 and endIdx or startIdx))
        if (direction == 1 and endIdx >= totalLength) or (direction == -1 and endIdx <= 1) then break end
        startIdx = endIdx + (direction == 1 and -1 or 1) * math.max(1, math.floor(currentSelectionLength * overlapFraction + 0.5))
    end
    return selectionStartArr, selectionEndArr
end

--------------------------------------------------------Number of rebuilds
--Calculate the number of rebuilds for each protein depending on Activity map
-- Function to calculate total activity for a selection
function calculateSelectionActivity(activityMap, startIdx, endIdx)
    local totalActivity = 0
    for i = startIdx, endIdx do
        totalActivity = totalActivity + activityMap[i]
    end
    return totalActivity
end

-- Function to distribute rebuilds based on activity and number of elements
function DistributeRebuildsByActivity(activityMap, numOfRebuilds)
    local selectionActivities = {}
    local selectionLengths = {}
    local totalActivity = 0
    local totalLength = 0

    -- Calculate activity and length for each selection
    for i = 1, #selectionStartArr do
        local activity = calculateSelectionActivity(activityMap, selectionStartArr[i], selectionEndArr[i])
        local length = selectionEndArr[i] - selectionStartArr[i] + 1
        table.insert(selectionActivities, activity)
        table.insert(selectionLengths, length)
        totalActivity = totalActivity + activity

        totalLength = totalLength + length
    end

    -- Calculate number of rebuilds for each selection based on activity
    local activityBasedRebuilds = {}
    local distributedActivityRebuilds = 0
    local numOfActivityRebuilds = math.floor(numOfRebuilds * 2 / 3)  --distribute just 67% of rebuilds


    for i = 1, #selectionStartArr do
        local proportion = selectionActivities[i] / totalActivity
        local rebuilds = math.floor(proportion * numOfActivityRebuilds)
        table.insert(activityBasedRebuilds, rebuilds)
        distributedActivityRebuilds = distributedActivityRebuilds + rebuilds
		--print ("+++++, i, rebuilds, proportion", i, rebuilds, proportion, numOfActivityRebuilds, numOfRebuilds)
    end


    -- Calculate number of rebuilds for each selection based on length
    local lengthBasedRebuilds = {}
    local distributedLengthRebuilds = 0
    local numOfLengthRebuilds = numOfRebuilds - numOfActivityRebuilds


    for i = 1, #selectionStartArr do
        local proportion = selectionLengths[i] / totalLength
        local rebuilds = math.floor(proportion * numOfLengthRebuilds)
        table.insert(lengthBasedRebuilds, rebuilds)
        distributedLengthRebuilds = distributedLengthRebuilds + rebuilds
    end

    -- Combine the two distributions
    local rebuildsDistribution = {}
    for i = 1, #selectionStartArr do
        rebuildsDistribution[i] = activityBasedRebuilds[i] + lengthBasedRebuilds[i]
    end

    -- Adjust for any rounding errors to ensure the total number of rebuilds equals numOfRebuilds
    local distributedRebuilds = distributedActivityRebuilds + distributedLengthRebuilds
    local remainingRebuilds = numOfRebuilds - distributedRebuilds
    --print(numOfRebuilds .. ":", "Remaining Rebuilds:", remainingRebuilds .. "-" .. distributedRebuilds .. "=" .. distributedRebuilds)


    while remainingRebuilds > 0 do
        local randomIndex = math.random(1, #rebuildsDistribution)
        rebuildsDistribution[randomIndex] = rebuildsDistribution[randomIndex] + 1
        remainingRebuilds = remainingRebuilds - 1
    end

    return rebuildsDistribution
end


--------------------------------------------------------service funtions

-- Define a function to calculate energy score per number of segments for the selection
function CalculateSelectionEnergyPerSegment(start, finish)
	local segmentCount = finish - start + 1
	local energyScore = 0
	for i = start, finish do
		--energyScore = energyScore + current.GetSegmentEnergyScore(i)
		energyScore = energyScore + GetSegmentBBScore(i)
	end
	return energyScore / segmentCount
end

-- Define a function to sort selections based on energy score per number of segments
function SortSelections(selectionStartArr, selectionEndArr)
	-- Create a table to store selection indices, their corresponding energy scores per segment, start index, and end index
	local sortedSelections = {}

	-- Calculate and store energy scores per segment for each selection
	for i = 1, #selectionStartArr do
		local start = selectionStartArr[i]
		local finish = selectionEndArr[i]
		local energyScorePerSegment = CalculateSelectionEnergyPerSegment(start, finish)
		table.insert(sortedSelections, {index = i, start = start, finish = finish, scorePerSegment = energyScorePerSegment})
	end

	-- Sort the sortedSelections table based on energy scores per segment in ascending/descending order
	if worstFirst then
		table.sort(sortedSelections, function(a, b) return a.scorePerSegment < b.scorePerSegment end) --ascending
	else
		table.sort(sortedSelections, function(a, b) return a.scorePerSegment > b.scorePerSegment end) --descending 
	end

	-- Print the start index, end index, energy score per segment of each selection, and the sorted start and end indices
	for _, selection in ipairs(sortedSelections) do
		if reportLevel>3 then print("Start index:", selection.start, "- End index:", selection.finish, "- Energy", selection.scorePerSegment) end
	end

	-- Return the sorted arrays
	local sortedStartArr = {}
	local sortedEndArr = {}
	for _, selection in ipairs(sortedSelections) do
		table.insert(sortedStartArr, selection.start)
		table.insert(sortedEndArr, selection.finish)
	end

	return sortedStartArr, sortedEndArr
end

--returns the id of the top X entries with the highest scoreBB values from the remixBBScores array:
function returnXinOrder(remixBBScores, X)
	SortByBackbone()
	-- Sort the array based on scoreBB in descending order
	--table.sort(remixBBScores, function(a, b) return a.scoreBB > b.scoreBB end)

	local topIds = {}
	for i = 1, math.min(X, #remixBBScores) do
		table.insert(topIds, remixBBScores[i].id)
	end
	if reportLevel>3 then print("topIds ".. #topIds..":", table.concat(topIds, " ")) end 
	return topIds
end

----------------------------------------------------------
  function Cleanup(err)
	print ("Cleanup")
    currentScore=ScoreReturn()
	behavior.SetClashImportance(endCI)
    undo.SetUndo(true) 

	save.Quickload(1)
	print ("Slot 1 score", ScoreReturn(),"bb",ScoreBBReturn())
	if lastSlots~=nil and false then --disabled
		for i=1, #lastSlots do
			save.Quickload(lastSlots[i])
			print ("Slot",lastSlots[i],"score", ScoreReturn(),"bb",ScoreBBReturn())
		end
	end

	save.Quickload(100)
	print ("Slot 100 score", ScoreReturn(),"bb",ScoreBBReturn())
    save.Quickload(98)
	print ("Slot 98 score", ScoreReturn(),"bb",ScoreBBReturn())

	if reportLevel>1 then  print ("Total Rebuilds Tested:", totalRebuildsDone) end

    --if convertLoop then save.LoadSecondaryStructure() end
	print (err)
  end

 -----------------------------------------------------------------------------------------------FUZE----------------------------------------------------------------------------------------------
function makeShake()
    behavior.SetClashImportance(1)
    structure.ShakeSidechainsAll (1)
end
function TinyFuze()
    SelectionSphere()
    behavior.SetClashImportance(0.05)
    structure.ShakeSidechainsSelected(1)
    behavior.SetClashImportance(1*CI)
    SetSelection()
end

---

-- Function to parse the input string into a 2D array
function parseInput(input)
    local result = {}
    for quadruplet in input:gmatch("{([^}]+)}") do
        local values = {}
        for value in quadruplet:gmatch("[^,%s]+") do
            table.insert(values, tonumber(value))
        end
        table.insert(result, {
            clashImportance = values[1],
            shakeIter = values[2],
            clashImportance2= values[3],
            wiggleIter = values[4]
        })
    end
    return result
end
-- Universal function that runs the logic on parsed input
function FuzeLocal(inputString)
	local startTime3 = os.clock()
	save.Quicksave(99) 

    local fuzeConfig = parseInput(inputString)
    local score = ScoreReturn()

    for idx, triplet in ipairs(fuzeConfig) do
        behavior.SetClashImportance(triplet.clashImportance)
        if triplet.shakeIter and triplet.shakeIter > 0 then
            structure.ShakeSidechainsAll(triplet.shakeIter)
        elseif triplet.shakeIter and triplet.shakeIter < 0 then
            structure.LocalWiggleAll(-triplet.shakeIter, false, true)
        end
        behavior.SetClashImportance(triplet.clashImportance2 )
        if triplet.wiggleIter and triplet.wiggleIter > 0 then
            structure.LocalWiggleAll(triplet.wiggleIter)
        elseif triplet.wiggleIter and triplet.wiggleIter < 0 then
            structure.LocalWiggleAll(-triplet.wiggleIter, true, false)
        end

		if triplet.clashImportance2==1 then 
			if  ScoreReturn() > score  then
				score =  ScoreReturn()
				save.Quicksave(99)
			else
				save.Quickload(99)
			end
		end

    end

	local endTime3 = os.clock()
	executionTime3 = executionTime3 + endTime3 - startTime3
    return score -- Return the score
end
---


function FuzeLocalGlobal1delete()
    behavior.SetClashImportance(0.05)
    structure.ShakeSidechainsAll(1) 
    structure.LocalWiggleAll(5)	
    
    behavior.SetClashImportance(1) 
	structure.ShakeSidechainsAll(2)
    undo.SetUndo(true) 
    structure.LocalWiggleAll(20)
    undo.SetUndo(false) 
end
function FuzeLocalGlobal2()
    behavior.SetClashImportance(0.05)
    structure.ShakeSidechainsAll(1) 
    structure.WiggleAll(5)	
    
    behavior.SetClashImportance(1) 
	structure.ShakeSidechainsAll(2)
    undo.SetUndo(true) 
    structure.WiggleAll(20)
    undo.SetUndo(false) 
end
function Fuze_low() --not used
	currentScore1=ScoreReturn()

	behavior.SetClashImportance(0.05)
	structure.ShakeSidechainsAll (1)
	structure.WiggleAll (5)

	behavior.SetClashImportance(0.05)
	structure.WiggleAll (3)
	behavior.SetClashImportance(1*CI)
	undo.SetUndo(true)
	structure.WiggleAll (20)
	undo.SetUndo(false) 
	if reportLevel > 3 then print ("fuze finished at", ScoreReturn()) end
	if convertLoop then save.LoadSecondaryStructure() end
end

function Fuze2()
	currentScore1=ScoreReturn()

	if bestSlot==0 then 
		bestSlot=99 
		print ("++++++++++++++++NO BESTSLOT VALUE!!!!!!!!!!!!!")
	end

	behavior.SetClashImportance(1)
	structure.ShakeSidechainsAll (1)
	behavior.SetClashImportance(0.05)
	structure.WiggleAll (10)

	behavior.SetClashImportance(0.05)
	structure.ShakeSidechainsAll (2)
	behavior.SetClashImportance(1)
	structure.WiggleAll (20)

	if ScoreReturn() > currentScore1 then 
	  save.Quicksave(bestSlot) 
	else save.Quickload(bestSlot) end

	if reportLevel > 3 then print ("fuze finished at", ScoreReturn()) end
	if convertLoop then save.LoadSecondaryStructure() end
end

function Fuze()
	currentScore1=ScoreReturn()
	
	save.Quicksave(99) 
	behavior.SetClashImportance(0.05)
	structure.ShakeSidechainsAll (1)
	structure.WiggleAll (2)

	behavior.SetClashImportance(0.25)
	structure.WiggleAll (2, 0, 1)
	structure.WiggleAll (2)

	behavior.SetClashImportance(1)
	structure.ShakeSidechainsAll (1)
	undo.SetUndo(true) 
	structure.WiggleAll (5, 1, 0)
	undo.SetUndo(false) 
	
	if ScoreReturn() > currentScore1 then 
	  save.Quicksave(99) 
	  currentScore1=ScoreReturn()
	end

	behavior.SetClashImportance(0.25)
	structure.ShakeSidechainsAll (1)
	structure.WiggleAll (2)

	behavior.SetClashImportance(1)
	undo.SetUndo(true)
	structure.WiggleAll (20)
	undo.SetUndo(false) 

	if ScoreReturn() > currentScore1 then 
	  save.Quicksave(99) 
	  currentScore1=ScoreReturn()
	end

	if reportLevel > 2 then print ("fuze finished at", ScoreReturn(), "bb", ScoreBBReturn()) end
	if convertLoop then save.LoadSecondaryStructure() end
end

---------------------------------------------------------------General service functions
function roundX(x, digits)
  digits = digits or 1
  local multiplier = 10 ^ digits
  return math.floor(x * multiplier) / multiplier
end

--[[ function returning the max between two numbers --]]
function max(num1, num2)
   if (num1 > num2) then  result = num1
   else result = num2    end
   return result
end

--	BackBone-Score is just the general score without clashing. Clashing is usefull to ignore when there is need to rank a lot of the Rebuild solutions very fast without the Fuze.
function ScoreBBReturn()
    x = 0
    for i=1, proteinLength do
      x = x + current.GetSegmentEnergySubscore(i, "Clashing")
    end
    x = current.GetEnergyScore() - x
	return x-x%1 --no decimals
end

-- Create array of the Scores for every Subscore of the puzzle
function GetSolutionSubscores(SolutionID)
    local scoreParts = puzzle.GetPuzzleSubscoreNames()
    local solutionSubscores = {}

	-- Iterate over each score part
	for _, scorePart in ipairs(scoreParts) do
		local currentSubscore = 0
		-- Iterate over each segment and calculate subscore
		for segmentIndex = 1, proteinLength do
			currentSubscore = currentSubscore + current.GetSegmentEnergySubscore(segmentIndex, scorePart)
		end
		solutionSubscores[scorePart] = currentSubscore
	end
    return solutionSubscores
end

function GetSegmentBBScore(i)
	return current.GetSegmentEnergyScore(i) - current.GetSegmentEnergySubscore(i, "Clashing")
end

-------------------------------------------------------------------------Selection functions

--selections
function FindSelection()
	selNum=0
	for k=1, proteinLength do
	  if selection.IsSelected(k) then 
		--find for selection start
		if k==1 then 
		  selNum=selNum+1
		  selectionStartArr[selNum]=k
		else  
		  if not selection.IsSelected(k-1) then 
			selNum=selNum+1
			selectionStartArr[selNum]=k 
		  end
		end
		
		--find the selection end
		if k==proteinLength then
		  selectionEndArr[selNum]=k
		else
		  if not selection.IsSelected(k+1) then selectionEndArr[selNum]=k end
		end
	   
	  end
	end

	--if no selection, select residues 3-6
	if selNum==0 then
	  selectionStartArr[1]=math.min(5,proteinLength)
	  selectionEndArr[1]=math.min(13,proteinLength)
	end
end

function SetSelection()
	selection.DeselectAll()
	for k=1, proteinLength do
		if (k>=selectionStart) and (k<=selectionEnd) then
			selection.Select (k)
		end
	end
end

function SelectLoops()
	selection.DeselectAll()
	looplength=0

	--select all loops with length >= 3
	for k=1, proteinLength do
	  --print (structure.GetSecondaryStructure(k))
	  if structure.GetSecondaryStructure(k) == 'L' then 
		looplength=looplength+1
		selection.Select (k) 
	  else
		if (looplength>0)  and (looplength<3) then  --deselect if loop is too short
		  for j=1, looplength do
			selection.Deselect (k-j)
		  end
		end
		looplength=0
	  end
	  
	  if (k==proteinLength) and (looplength==1) then selection.Deselect (k) end
	end
end
  
function SetAllSelections()
	selection.DeselectAll()
	for j=1, selNum do
	  selectionStart=selectionStartArr[j]
	  selectionEnd=selectionEndArr[j]
	  if (selectionEnd-selectionStart < 0) then
		temp = selectionStart
		selectionStart = selectionEnd
		selectionEnd = temp
	  end
	  if (reportLevel>3) then print ("Setting selection"..j.."/"..selNum..": ", selectionStart.."-"..selectionEnd) end
	  for k=1, proteinLength do
		if (k>=selectionStart) and (k<=selectionEnd) then
		  selection.Select (k)
		end
	  end
	end
end

-- Select everything in radius sphereRadius(=8) near any selected segment.
function SelectionSphere()
	--dump selection to array
	selectedSegs={}
	for k=1, proteinLength do
	  if selection.IsSelected(k) then  
		selectedSegs[k]=1
	  else
		selectedSegs[k]=0
	  end
	end

	for k=1, proteinLength do
	  for j=1, proteinLength do
		dist_str = structure.GetDistance(k, j)
		if (selectedSegs[j] == 1) and (dist_str < sphereRadius) then
		   selection.Select(k)
		end
	  end
	end
end

-- Function to print selections with the number of rebuilds for each
function printSelections()
	if (reportLevel>3) then print("Creating "..#rebuildsDistribution.." selections") end
    for i = 1, #selectionStartArr do
        local rebuilds = rebuildsDistribution[i]
        if (reportLevel>2) then print("Selection " .. i .. ": " .. selectionStartArr[i].." - "..selectionEndArr[i]," Number of Rebuilds: " .. rebuilds) end
    end
end



----------------------------------------------------------------------------------------Inheritance statistics system---------------------------------------------------------------------------------
history = {}
historySize = 0
uniqueID = 1
historyIDbySlot = {}
historyIDbySlot[1] = 1
lastBestSlot = 0
newMask = {}
for i = 1, proteinLength do table.insert(newMask, 0) end
table.insert(history, {id = 1, slot = 1, parent = nil, offspring = {}, scoreBB = ScoreBBReturn(), score = ScoreReturn(), changeMask = newMask, generation = 0}) 
--history[currentID] = {id = currentID, slot = slot, parent = parentID, offspring = {}, scoreBB = scoreBB, changeMask = mask or {}}

function addHistoryEntry(slot, parentSlot, scoreBB, mask, generation)
	--if we fill the array after the re	set then use last bestslot as a parent
	if parentSlot == 1 and lastBestSlot ~= 0 then
		parentSlot = lastBestSlot
	end
	uniqueID = uniqueID + 1
	historySize = historySize + 1
	local currentID = uniqueID
	parentID = historyIDbySlot[parentSlot]
	historyIDbySlot[slot] = currentID
	if not history[currentID] then
		history[currentID] = {id = currentID, slot = slot, parent = parentID, offspring = {}, scoreBB = scoreBB, score = 0, changeMask = mask, generation = generation}
    end
	if parentID then
		if parentID == currentID then
			if reportLevel > 2 then print("Error: A slot cannot be its own parent. Skipping invalid parent assignment.", parentSlot, parentID) end
			history[currentID].parent = nil
		else
			if not history[parentID] then
				history[parentID] = {id = parentID, slot = 0, parent = nil, offspring = {currentID}, scoreBB = nil, score = 0}
			end
			table.insert(history[parentID].offspring, currentID)
		end
	end
	return currentID
end

function updateHistoryArray (slot, currentScore, currentBBScore)
	idx = historyIDbySlot[slot]
	--history[idx].scoreBB = currentBBScore
	history[idx].score = currentScore
end
	
function resetInheritanceTable(slot)
	--we save bestSlot on reset to use it in addHistoryEntry() 
	--when filling table with parentSlot==1 (to use ==currentBestSlot instead to save the inheritance)
	lastBestSlot = slot
	local currentID = historyIDbySlot[slot]
	--remove everything except current and 1st ID
	for idx, info in pairs(history) do
		if idx~= 1 then
			history[idx].slot = 0
		end
	end
end
function removeSlotStats(slot)
	local id = historyIDbySlot[slot]
	if history[id] then
		history[id].slot = 0
	end
end
function removeHistoryByID(id)
	--first, remove offspring from the parentID record
	parentID = history[id].parent
	local offspring = history[parentID].offspring
    if #offspring > 0 then
        table.remove(offspring, #offspring)
    end
   historySize = historySize -1
   history[id] = nil
end
---------------------------------------Calculate Rarity
cumulativeChangeMasks = {} --store mask for each id to optimize speed
relationScores = {} -- store relation scores between pairs of elements
function combineMasks(mask1, mask2)
    local combined = {}
    for i = 1, #mask1 do
        combined[i] = math.max(mask1[i], mask2[i])
    end
    return combined
end
function calculateRelation(mask1, mask2)
    local totalChanges = #mask1
    local diffCount = 0
    for i = 1, totalChanges do
        if mask1[i] ~= mask2[i] then
            diffCount = diffCount + 1
        end
    end
    return 1 - (diffCount / totalChanges)
end
local function getCumulativeChangeMask(id)
    if cumulativeChangeMasks[id] then
        return cumulativeChangeMasks[id]
    end
    local changeMask = history[id].changeMask or {}
    local parentId = history[id].parent
    while parentId do
        local parentChangeMask = history[parentId].changeMask or {}
        changeMask = combineMasks(changeMask, parentChangeMask)
        parentId = history[parentId].parent
    end
    cumulativeChangeMasks[id] = changeMask
    return changeMask
end
local function getRelationScore(idA, idB)
    local key = idA < idB and idA .. "-" .. idB or idB .. "-" .. idA
    if relationScores[key] then
        return relationScores[key]
    end
    local maskA = getCumulativeChangeMask(idA)
    local maskB = getCumulativeChangeMask(idB)
    local relationScore = calculateRelation(maskA, maskB)
    relationScores[key] = relationScore
    return relationScore
end

function calculateRarityScore(id)
    local totalRelation = 0
    local count = 0
    for otherID, _ in pairs(history) do
		--print (id, otherID, history[otherID].slot)
        if (id ~= otherID) and (history[otherID].slot ~= 0) then
            local relationScore = getRelationScore(id, otherID)
            totalRelation = totalRelation + relationScore
            count = count + 1
        end
    end
    if count == 0 then  		return 1   end
    if totalRelation == 0 then 	return 1   end
	local rarityScore = count / totalRelation - 1
    if reportLevel > 3 then print("ID:", id, "Total Relation:", totalRelation, "Count:", count, "Average Relation:", averageRelation, "Normalized Relation:", normalizedRelation) end
    return rarityScore
end

----------------------------------------------------------
local function getAncestors(id)
    local ancestors = {}
    local visited = {} -- Track visited nodes to detect cycles
    while id do
        if visited[id] then
            if reportLevel > 2 then print("Cycle detected! Stopping ancestor retrieval to avoid infinite loop.") end
            break
        end
        table.insert(ancestors, id)
        visited[id] = true
        id = history[id] and history[id].parent or nil
        if reportLevel > 4 then print("Retrieving ancestor, current id:", id) end
    end
    return ancestors
end
function countDescendants(id)
    local function countRecursive(currentID)
--print ("currentID",currentID, id)
        if not history[currentID] then return 0 end -- Add this check to prevent crashes
        local count = 0
        if (#history[currentID].offspring > 0) then
			for _, offspringID in ipairs(history[currentID].offspring or {}) do
				if history[offspringID] then  -- Check if offspringID exists in history
					if history[offspringID].slot ~= 0 then  -- Do not count removed records
						count = count + 1
					end
					-- Continue recursion even if current isn't active, to find active descendants
					count = count + countRecursive(offspringID)
				else
					--print("!") UNCOMMENT TO FIND THE BUG!!!!!!!!!!!!!!!!!!!!!!!!!!!!
				end
			end
		end
        return count
    end
    return countRecursive(id)
end
function calculateDescendantPercentage(id)	
    local totalDescendants = countDescendants(id)
    local totalSlots = 0 
    for idx, _ in pairs(history) do
		if history[idx].slot ~= 0	then
			totalSlots = totalSlots + 1
		end
    end
	if totalSlots == 0 then totalSlots = 1 end -- prevent div zero
    return (totalDescendants / totalSlots) * 100
end

local function printHistoryTree(id, level)
    local info = history[id]
    --print (id, level, info.slot)
    if info then
        local indent = string.rep("\t", level)
		if info.slot ~= 0 then
			local parentID = info.parent or "None"
			local scoreBB = info.scoreBB

			--local rarity = calculateBreedRarityScore(id)
			local rarity = calculateRarityScore(id)
			remixBBID = getRemixBBIDBySlot (info.slot) 
			if remixBBID == 0 then rarity, rank, descriminator = 0, "", 0
			else 
				--rarity =  remixBBScores[remixBBID].rarity
				rank = remixBBID
				descriminator = remixBBScores[remixBBID].bbDescriminator
				descriminator =  (descriminator>=0) and string.format("+%.0f",descriminator) or string.format("%.0f",descriminator)
				if rank < 10 then rank = tostring (rank).."++"
				else rank = tostring (rank) end
			end

			local descendantPercentage = calculateDescendantPercentage(id)
			local descendantText = ""
			if descendantPercentage > 0 then descendantText = string.format(", Descendant%%: %.2f", descendantPercentage) end
			print(string.format("%s%sSlot %2d: ID: %3d, ParentID: %3s, ScoreBB: %.0f  %s, score: %.0f, Gen: %3d Rarity: %.2f%s", rank, indent, info.slot, id, parentID, scoreBB, descriminator, info.score, info.generation, rarity, descendantText ))

			-- Sort offspring by scoreBB before printing them
			table.sort(info.offspring, function(a, b)
				local scoreA = history[a] and history[a].scoreBB or 0
				local scoreB = history[b] and history[b].scoreBB or 0
				return scoreA > scoreB
			end)

			for _, offspringID in ipairs(info.offspring) do
				printHistoryTree(offspringID, level + 1)
			end
		else
			--if it has some active offsprings, draw an empty line instead and continue recursion
			if #info.offspring > 0 then
				print(indent..string.rep("- ", 15))
				for _, offspringID in ipairs(info.offspring) do
					printHistoryTree(offspringID, level + 1)
				end
			end
		end	--if info.slot
    end --if info 
end

-- Clean the history table from the records with no active offsprings
function cleanHistory()
	local toRemove = {}
	for id, info in pairs(history) do
		if (#info.offspring > 0) then
			local totalDescendants = countDescendants(id)
			if totalDescendants == 0 then
				if  (info.slot ==0) then
					table.insert(toRemove, id)
				else
					history[id].offspring = {}  
				end
			end
		end
	end
	-- Remove collected IDs
	for _, id in ipairs(toRemove) do
		history[id] = nil
	   historySize = historySize -1
    end
end

counterx = 0
-- Prints the entire history in a tree structure
function printHistory()
	cleanHistory()
	print ("Print history tree ("..historySize..") : ")
	--search for the first record to start recursion
    for id, info in pairs(history) do
        if not info.parent then
            printHistoryTree(id, 0)
        end
    end
end

-------------------------------------------------------------


xpcall ( main , Cleanup )