Code
-- jeff101 began DistMap0 12/3/15 354pm
-- adapted DistMap0c.txt 12/4/15 1220pm noon code to make DistMap1
-- last updated 1/24/16 312am
--
-- version 1.1 LociOiling 2016/02/14
-- + double sample depth to two hex digits (8 bits) per sample
--
-- global variables
recipename='DistMap1'
function Main()
local dist,dd,dcuts,j,j1,j2,minj1,minj2,maxj1,maxj2,ask,runflag,isok
local chlist,nshades,hnshades,chway,revflagi,revflag,chways,spflag,lbls
local sss,aas,dists,dds,mindist,maxdist,tmpstr,epsway,epsdim,makecsv,maketot
local csvspacing,dtot,dtotcalc,davg,drms,drmsd,dmed,dtopq,dbotq,dlist
local justansi,ansimax,nlists,didit,skip5num,pflag,sortway
local tot=structure.GetCount()
local htot=math.floor(tot/2) -- rounds half the total down to nearest integer
local nbands=band.GetCount()
print(' ')
print(recipename..': Puzzle '..puzzle.GetName())
print(' with Puzzle ID '..puzzle.GetPuzzleID()..' on '..ui.GetPlatform()..' at '..os.date())
print(' has score '..trunc3(current.GetScore())..', '..nbands..' bands, and '..tot..' residues.')
print(' ')
chlist=makechlist() -- chlist should contain many lists chlist[j] containing 5 or 10 elements each
nlists=#chlist -- # of character sets available
ansimax=nlists-2 -- lists 1 to ansimax are ANSI-only,
-- lists ansimax+1 to nlists contain non-ANSI chars
--
-- below are some default values
-- first few hold for most charts
chway=1 -- which set of characters to use for shading the charts
revflagi=false -- false=regular order, true=reverse order for shades
justansi=true -- true=just use ansi chars, false=allow non-ansi chars
spflag=0 -- 0=no spaces, 1=single space, 2=commas
lbls=2 -- 0=all labels off, 1=left labels off, top/bottom labels on, 2=all labels on
maketot=false -- false=no full-size chart, true=makes full-size chart
-- next few are for csv charts
makecsv=3 -- 0=none, 1=1 csv chart, 2=2 csv charts, 3=3 csv charts
csvspacing=2 -- 1=single space, 2=commas
-- next few are for eps plot
epsway=4 -- 0=no eps plot, 1=5-shade eps plot, 2-4=10-shade eps plot (4>3>2 for short-distance contrast)
epsdim=tot -- size of eps plot to make
-- next few are general ones
skip5num=0 -- 0=does all 5-number-summaries, 1=skips them for groups >= 500 elements, 2=skips for all groups
sortway=2 -- 1=non-table.sort() way, 2=table.sort() way
--
showchlist(chlist,nlists,ansimax,1) -- show ANSI-only sets in chlist
didit=0
--
-- below shows input menus until the inputs make sense
isok=0 -- will be 1 when the inputs make sense
while isok==0 do
ask=dialog.CreateDialog(recipename)
ask.LabelA=dialog.AddLabel('------------- VALUES FOR MOST CHARTS --------------')
if justansi then
ask.chway=dialog.AddSlider("chway: ",chway,1,ansimax,0)
else
ask.chway=dialog.AddSlider("chway: ",chway,1,nlists,0)
end
tmpstr='chway picks characters for distance ranges in charts.\nRecipe Output and scriptlog.default.xml give details.'
ask.LabelB=dialog.AddLabel(tmpstr)
ask.revflagi=dialog.AddCheckbox('Reverse the order of characters in most charts.',revflagi)
ask.justansi=dialog.AddCheckbox('Use only ANSI chars (makes smaller output files).',justansi)
ask.spflag=dialog.AddSlider('chart spacing: ',spflag,0,2,0)
ask.LabelC=dialog.AddLabel('spacing: 0=no spaces, 1=single space, 2=comma.')
ask.lbls=dialog.AddSlider('chart labels: ',lbls,0,2,0)
ask.LabelD=dialog.AddLabel('labels: 0=all off, 1=left off, top/bot on, 2=all on.')
ask.maketot=dialog.AddCheckbox('Include a full-size chart.',maketot)
ask.LabelE=dialog.AddLabel('--------- VALUES FOR EXCEL/CSV CHARTS ---------')
ask.makecsv=dialog.AddSlider('makecsv: ',makecsv,0,3,0)
ask.LabelF=dialog.AddLabel('makecsv: # of Excel/csv charts to make.\n Each successive chart is larger.')
ask.csvspacing=dialog.AddSlider('csv spacing: ',csvspacing,1,2,0)
ask.LabelG=dialog.AddLabel('spacing: 1=single space, 2=comma.')
ask.LabelH=dialog.AddLabel('----------------- VALUES FOR EPS PLOT -----------------')
ask.epsway=dialog.AddSlider('epsway: ',epsway,0,4,0)
ask.LabelI=dialog.AddLabel('epsway: 0=no plot, 1=5 shades, 2-4=10 shades.\n epsway 4>3>2 for short-distance contrast.')
ask.epsdim=dialog.AddSlider('epsdim: ',epsdim,2,tot,0)
ask.LabelJ=dialog.AddLabel('epsdim: # of AA\'s to show in an eps plot.')
ask.LabelK=dialog.AddLabel('------------------ GENERAL CONTROLS -------------------')
-- ask.skip5num=dialog.AddSlider('skip5num: ',skip5num,0,2,0)
-- ask.LabelL=dialog.AddLabel('skip5num: skip 5-#-summaries (median, top & bot 1/4)\n 0=never, 1=for sets>=500 elements, 2=always.')
ask.sortway=dialog.AddSlider('sortway: ',sortway,1,2,0)
ask.LabelM=dialog.AddLabel('sortway: 1=slow way, 2=table.sort() way.')
ask.run=dialog.AddButton('Run',1)
ask.quit=dialog.AddButton('Quit',0)
runflag=dialog.Show(ask)
if runflag==0 then -- quit
isok=1
print('Quitting now.')
else
isok=1
-- values for most charts
chway=ask.chway.value+0
revflagi=ask.revflagi.value
if revflagi then
revflag=1
else
revflag=0
end -- if revflagi
if justansi~=ask.justansi.value then -- the ANSI setting has changed
if justansi then -- new setting allows non-ANSI chars
if didit==0 then
showchlist(chlist,nlists,ansimax,2) -- show the non-ANSI character sets in chlist
didit=1
chway=nlists -- reset the default value to a non-ANSI list
end -- if didit
isok=0 -- make sure to repeat the input menu
end -- if justansi
end -- if justansi
justansi=ask.justansi.value
if justansi then
if chway>ansimax then
print('ERROR: chway='..chway..' contains non-ANSI characters.')
print(' Please pick another character set or allow non-ANSI characters.')
print(' ')
chway=1 -- reset the default value to an ANSI character set
isok=0
end -- if chway
end -- if justansi
spflag=ask.spflag.value+0
lbls=ask.lbls.value+0
maketot=ask.maketot.value
-- csv values
makecsv=ask.makecsv.value+0
csvspacing=ask.csvspacing.value+0
-- eps values
epsway=ask.epsway.value+0
epsdim=ask.epsdim.value+0
-- general values
-- skip5num=ask.skip5num.value+0
sortway=ask.sortway.value+0
end -- if runflag
end -- if isok
--
if runflag==1 then
tmpstr=''
nshades=#chlist[chway]
for j=1,nshades do
tmpstr=(tmpstr..chlist[chway][j])
end -- for j
--
-- list the input values
print('Got the following inputs by '..os.date()..' for '..tot..' residues:')
print(' chway='..chway..'=\''..tmpstr..'\'='..nshades..' shades, reverse order='..tostring(revflagi)..'='..revflag..',')
print(' use only ANSI characters='..tostring(justansi)..', spacing='..spflag..',')
print(' labels='..lbls..', make full-size chart='..tostring(maketot)..',')
print(' makecsv='..makecsv..', csvspacing='..csvspacing..', epsway='..epsway..', epsdim='..epsdim..',')
print(' skip5num='..skip5num..', and sortway='..sortway..'.')
print(' ')
--
-- next part calculates Ca-Ca distance map in dists and dds
--
print('Gathering data and doing calculations...')
sss={}
aas={}
dists={}
dds={}
mindist= -1 -- will be min nonzero Ca-Ca distance found
maxdist= -1 -- will be max nonzero Ca-Ca distance found
dcuts={2,4,6,8,10,15,20,30,50} -- borders of distance ranges
-- want an odd # of elements in dcuts, preferably 9
for j1=1,tot do
dists[j1]={}
dds[j1]={}
sss[j1]=structure.GetSecondaryStructure(j1)
aas[j1]=structure.GetAminoAcid(j1)
end -- for j1
print('Got list of '..tot..' AA\'s and their secondary structures.')
dtot=0 -- will be total # of distances between 2 different Ca's
dlist={}
for j1=1,tot do
-- print('Getting distances from alpha-carbon '..j1..' of '..tot..' ('..os.date()..').')
for j2=j1,tot do -- j2>=j1 holds here
if j1==j2 then
dist=0+0
dd=getdd(dcuts,dist)
dists[j1][j2]=dist
dds[j1][j2]=dd
else -- j1~=j2, which means j2>j1
dtot=dtot+1
dist=getdist(j1,j2,0,0)
dlist[dtot]=dist
-- dist should now be distance between alpha-carbons j1 and j2
if dist>0 then
if maxdist== -1 or dist>maxdist then
maxdist=dist
maxj1=j1
maxj2=j2
end -- if maxdist
if mindist== -1 or dist<mindist then
mindist=dist
minj1=j1
minj2=j2
end -- if mindist
end -- if dist
dd=getdd(dcuts,dist)
dists[j1][j2]=dist
dists[j2][j1]=dist
dds[j1][j2]=dd
dds[j2][j1]=dd
end -- if j1==j2
end -- for j2
end -- for j1
print('Got '..dtot..' distances among '..tot..' alpha-carbons ('..os.date()..').')
--
-- below calcs some stats for all the above distances
dtotcalc=tot*(tot-1)/2 -- how many pairs of Ca's above give distances to average
if dtot~=dtotcalc then
print('ERROR: dtot='..dtot..' is not '..dtotcalc..' like it should be.')
end -- if dtot
--
print(' ')
test5num(skip5num,sortway) -- test get5num
--
-- below gets 5-number summary as in https://en.wikipedia.org/wiki/Five-number_summary
dtopq,dmed,dbotq,davg,drms,drmsd=get5num(dlist,dtot,maxdist,mindist,skip5num,sortway)
davgp=davg+drmsd -- avg + rms dev
davgm=davg-drmsd -- avg - rms dev
--
-- next part makes a few charts
--
-- below makes a few charts based on dists for use in Excel
if makecsv>0 then
csvcharts(sss,aas,dists,tot,htot,csvspacing,makecsv,skip5num,sortway)
end -- if makecsv
--
-- below makes an eps file
if epsway>0 then
epsmaker(chlist,epsway,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,davg,drmsd,tot,epsdim)
end -- if epsway
--
-- next part prints the distance map in dds
--
nshades=(#dcuts)+1 -- should be 10
hnshades=nshades/2 -- should be 5
print('Some of the alpha-carbon distance map(s) below use '..nshades..' shades:')
print(string.format(' 0 means distance is in range (-Inf,%2d).',dcuts[1]))
for j=1,nshades-2 do
print(string.format(' %d means distance is in range [%2d,%2d).',j,dcuts[j],dcuts[j+1]))
end -- for j
print(string.format(' %d means distance is in range [%2d,Inf).',nshades-1,dcuts[nshades-1]))
print('Some of the alpha-carbon distance map(s) below use '..hnshades..' shades:')
print(string.format(' 0 means distance is in range (-Inf,%2d).',dcuts[2]))
for j=1,hnshades-2 do
print(string.format(' %d means distance is in range [%2d,%2d).',j,dcuts[2*j],dcuts[2*(j+1)]))
end -- for j
print(string.format(' %d means distance is in range [%2d,Inf).',hnshades-1,dcuts[2*(hnshades-1)]))
print(string.format('All nonzero distances range from %.2f to %.2f',mindist,maxdist))
print(string.format(' for the Ca-Ca pairs %d-%d and %d-%d respectively.',minj1,minj2,maxj1,maxj2))
print(string.format('For these distances, rms value=%.2f, average=%.2f, rms deviation=%.2f,',drms,davg,drmsd))
print(string.format(' so max=%6.2f, avg+rmsdev=%6.2f, avg=%6.2f, avg-rmsdev=%6.2f, min=%6.2f,',maxdist,davgp,davg,davgm,mindist))
print(string.format(' and max=%6.2f, upper 1/4=%6.2f, med=%6.2f, lower 1/4=%6.2f, min=%6.2f.',maxdist,dtopq,dmed,dbotq,mindist))
-- the last line comes from https://en.wikipedia.org/wiki/Five-number_summary
print(' ')
if maketot then
makechart(chlist,chway,revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,tot,0)
print(' ')
end -- if maketot
if tot>74 then
makechart(chlist,chway,revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,74,0)
print(' ')
end -- if tot
if tot>34 then
makechart(chlist,chway,revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,34,0)
print(' ')
end -- if tot
if tot>18 then
makechart(chlist,chway,revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,18,0)
print(' ')
end -- if tot
if tot>9 then
makechart(chlist,chway,revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,9,0)
print(' ')
end -- if tot
if tot>=9 then -- below makes special-format charts for Recipe Output window
-- Recipe Output window's slider is hard to line up on these charts.
-- It automatically lines up with final chart if there is no blank line at the end.
-- Usually I put blank lines after each item I print, but here I put them before.
-- pflag is used below to make no blank line after final chart.
pflag=0 -- 0 means no printing yet, 1 means something has printed already.
--
-- below lists the all-ANSI character sets
-- chways={2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,19,20,21,22,1,9,10}
chways={1} -- use 1 to give output like in DistMap0
-- chways={1,9,10,24} have fixed-width characters and should look best in Recipe Output window
for chway=1,#chways do
if pflag==1 then print(' ') end -- if pflag
makechart(chlist,chways[chway],0,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,9,1)
pflag=1
end -- for chway
if not justansi then -- non-ANSI char sets are allowed
-- below lists all the non-ANSI character sets
-- chways={23,24} -- 23 looks odd in Recipe Output window
chways={24} -- 24 looks fine in Recipe Output window, so use it here
for chway=1,#chways do
for revflag=0,1 do -- do it normal order (for white bckgrnd windows), then in reverse order (for black bckgrnd windows)
if pflag==1 then print(' ') end -- if pflag
makechart(chlist,chways[chway],revflag,spflag,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,9,1)
pflag=1
end -- for revflag
end -- for chway
end -- if justansi
end -- if tot
end -- if runflag
end -- Main()
function makechlist()
local chlist={}
-- in below, most characters have same width in Courier font, diff widths in Arial.
-- all integers have same width on Recipe Output window (Arial not Courier).
-- all of the utf8(9608 9619 9618 9617) have same width on Recipe Output window,
-- but utf8(9615) (white box) has a slightly different width.
chlist[1]={'0','1','2','3','4','5','6','7','8','9'} -- decimal
chlist[2]={'0','2','3','5','7','8','A','C','D','F'} -- hexadecimal for eps format
chlist[3]={'0','2','4','6','8','A','C','D','E','F'} -- hexadecimal for eps format
chlist[4]={'0','3','6','8','A','B','C','D','E','F'} -- hexadecimal for eps format
-- chlist[2] is evenly-spaced number-wise, but the low #'d shades are hard to tell apart
-- chlist[3-4] space the low #'d shades out more to raise contrast at short distances (4>3>2)
--
-- http://www.cs.umd.edu/Outreach/hsContest98/questions/node4.html gives
chlist[5]={'$','E','F','L','l','v','!',';',',','.'}
-- https://dboikliev.wordpress.com/2013/04/20/image-to-ascii-conversion/ gives
chlist[6]={'@','%','#','x','+','=',':','-','.',' '}
-- http://paulbourke.net/dataformats/asciiart/ gives
chlist[7]={'@','%','#','*','+','=','-',':','.',' '}
-- http://members.optusnet.com.au/astroblue/grey_scale.txt gives 12-char list
chlist[8]={'#','@','8','O','o','\"',';',',','.',' '} -- I omit 2 of them
--
chlist[9]={'0','1','2','3','4'} -- decimal
chlist[10]={'0','2','4','6','8'} -- decimal
chlist[11]={'0','4','8','B','F'} -- hexadecimal for eps format
-- next lines do order darkest to lightest
chlist[12]={'@','O','o','.',' '}
chlist[13]={'@','0','o','.',' '}
chlist[14]={'#','t','+','-',' '}
chlist[15]={'#','t','=','-',' '}
chlist[16]={'#','\"','~','-',' '} -- \" and utf8(34) are " character
chlist[17]={'#','\"','\'','`',' '} -- \' and utf8(39) are ' character
chlist[18]={'#','X','/',',',' '} -- utf8(96) is ` character
chlist[19]={'%','8','b','o',' '} -- \\ and utf8(92) are \ character
chlist[20]={'%','B','b','o',' '} -- utf8(47) is / character
chlist[21]={'%','8','b','c',' '}
chlist[22]={'%','B','b','c',' '}
chlist[23]={utf8(9608),utf8(9619),utf8(9618),utf8(9617),'_'}
chlist[24]={utf8(9608),utf8(9619),utf8(9618),utf8(9617),utf8(9472)}
--
-- in chlist[23], utf8(9615)=(white block) has odd width vs the rest
-- in both Arial and Courier fonts!
-- '_' at least has same width as #'s in both fonts
--
-- utf8(9472) and utf8(9552) are same width as shaded blocks in Arial & Courier fonts
-- utf8(9472) looks like -, utf8(9552) looks like =, and both are fairly white
-- they work better than the '_' character
--
-- http://www.w3schools.com/charsets/ref_utf_block.asp
-- http://www.sorenwinslow.com/AltCodesRun.asp and
-- https://en.wikipedia.org/wiki/Box-drawing_character give above codes
-- http://www.asciitable.com/ says utf8(176)-utf8(178) are shaded blocks too
--
-- typing raw shaded block characters above requires saving this file as unicode not ANSI.
-- saving the file as unicode takes about twice the space as ANSI (UTF-8 files = same size as ANSI ones).
-- will unicode recipe files put into the cookbook make all.macro much larger?
-- will they act strangely in the cookbook?
-- unicode outputs in scriptlog.default.xml files will also take twice the usual space.
--
return chlist
end -- makechlist()
-- chflag==1 prints ANSI-only char sets, chflag==2 prints sets w/non-ANSI chars
function showchlist(chlist,nlists,ansimax,chflag)
local j,k,str,chs,nchs
if chflag==1 then
print('Use chway to pick which character set (see chlist[chway] below) to use.')
print('Each character will represent a different distance range in the charts.')
print('Some character sets have 10 distance ranges while others have only 5.')
print('Some sets contain non-ANSI characters, which use more bytes than ANSI ones,')
print(' making scriptlog.default.xml larger for the same total # of characters.')
print(' ')
print('The following character sets use only ANSI characters:')
else -- if chflag==2 then
print('The following character sets include non-ANSI characters:')
end -- if chflag
-- all chlist[j]'s with j<=ansimax use ANSI characters only
for j=1,nlists do
if (j<=ansimax and chflag==1) or (j>ansimax and chflag==2) then
str=''
chs=chlist[j]
nchs=#chs
for k=1,nchs do
str=(str..chs[k])
end -- for k
print(string.format(' chlist[%2d][0-%d]=%2d shades=\'%s\'.',j,nchs-1,nchs,str)) -- \' prints a ' character
end -- if j
end -- for j
print(' ')
end -- showchlist()
-- code below came from
-- http://stackoverflow.com/questions/7983574/how-to-write-a-unicode-symbol-in-lua
function utf8(decimal)
if decimal<128 then return string.char(decimal) end
local bytemarkers = { {0x7FF,192}, {0xFFFF,224}, {0x1FFFFF,240} }
local charbytes = {}
local bytes,vals,b
for bytes,vals in ipairs(bytemarkers) do
if decimal<=vals[1] then
for b=bytes+1,2,-1 do
local mod = decimal%64
decimal = (decimal-mod)/64
charbytes[b] = string.char(128+mod)
end -- for b
charbytes[1] = string.char(vals[2]+decimal)
break
end -- if decimal
end -- for bytes
return table.concat(charbytes)
end -- utf8()
-- below finds the distance between residue r1 atom a1 and residue r2 atom a2
function getdist(r1,r2,a1,a2)
local tot,nuband
local dist= -1+0 -- default value
if a1==a2 and a1==0 then
dist=structure.GetDistance(r1,r2)
-- above finds distance between atom 0's, the alpha-carbons on each residue
else
-- below finds distance between general atoms on each residue
tot=band.GetCount()
nuband=band.AddBetweenSegments(r1,r2,a1,a2)
if nuband>tot then
dist=band.GetLength(nuband)
band.Delete(nuband)
else
-- 4/9/14 below happened for a1==a2 and a1==0 when r1 and r2 differed by 1
-- 5/26/14 below happened for a1==1 and a2==4 when r1 and r2 differed by 1 or less
print(string.format("ERROR: getdist(%d,%d,%d,%d) failed so returning -1.",r1,r2,a1,a2))
end -- if nuband
end -- if a1==a2
return dist
end -- getdist()
-- below returns an integer val from 0 to #dcuts
-- representing what range the distance dist lies in
-- val should be 0 if dist<dcuts[1] or (-Inf,dcuts[1])
-- val should be j if dcuts[j]<=dist<dcuts[j+1] or [dcuts[j],dcuts[j+1])
-- val should be #dcuts if dcuts[#dcuts]<=dist or [dcuts[#dcuts],Inf)
-- ideally #dcuts is 9 or less
function getdd(dcuts,dist)
local val=0+0
local j
for j=1,#dcuts do
if dist>=dcuts[j] then
val=j
end -- if dist
end -- for j
return val
end -- getdd()
-- below makes charts suitable for use in Excel
function csvcharts(sss,aas,dists,tot,htot,spflag,makecsv,skip5num,sortway)
local spc,j,k,m,js,jnk,gotit,tmpstr
local dist,davg,dmax,dmin,drms,drmsd,dstats,dlist,dtot,dtopq,dmed,dbotq
local maxj,maxk
--
if makecsv>0 then -- make at least one chart for Excel
--
-- usually spflag=0 no space, 1=single space, 2=commas
-- but here do 0,1=single space, 2=commas
--
if spflag==2 then
spc=','
else
spc=' '
end -- if spflag
--
-- below prepares and prints 1st chart, finding largest distance
-- remaining each time, then removes the 2 segments giving it,
-- like peeling an onion
gotit={}
for j=1,tot do
gotit[j]=0 -- 0 means not used yet, 1=used already
end -- for j
print('The next chart is a \'peel-the-onion\' chart. It is made as follows:')
print(' (1) Begin with a list of all alpha-carbons in the protein.')
print(' (2) Find which pair of alpha-carbons in the list is furthest apart')
print(' (find the outermost \'layer\' of the \'onion\').')
print(' (3) Report this pair and its members\' distance apart.')
print(' (4) Remove this pair from the list of alpha-carbons')
print(' (peel off the outermost \'layer\' of the \'onion\').')
print(' (5) Repeat steps (2)-(4) until no more alpha-carbon pairs')
print(' remain in the list.')
print(' (6) If there is one alpha-carbon left in the list,')
print(' report it as a pair with zero distance.')
print('In this chart, surface residues should have larger distances than core ones.')
print(' ')
print('********** the \'peel-the-onion\' chart below can be used in Excel ('..os.date()..') **********')
print(string.format('SS%sAA#%sAA%sSS%sAA#%sAA%s dist',spc,spc,spc,spc,spc,spc))
for m=1,htot do
dmax= -1
maxj=0
maxk=0
for j=1,tot-1 do
if gotit[j]==0 then
for k=j+1,tot do -- here k>j always holds
if gotit[k]==0 then
if dists[j][k]>dmax then
dmax=dists[j][k]
maxj=j
maxk=k
end -- if dists
end -- if gotit
end -- for k
end -- if gotit
end -- for j
j=maxj
k=maxk
gotit[j]=1
gotit[k]=1
print(string.format(' %s%s%3d%s %s%s %s%s%3d%s %s%s%6.2f',sss[j],spc,j,spc,aas[j],spc,sss[k],spc,k,spc,aas[k],spc,dists[j][k]))
end -- for m
if 2*htot<tot then -- tot is odd, so only 1 segment should remain w/gotit[j]==0
for j=1,tot do
if gotit[j]==0 then
k=j
gotit[j]=1
print(string.format(' %s%s%3d%s %s%s %s%s%3d%s %s%s%6.2f',sss[j],spc,j,spc,aas[j],spc,sss[k],spc,k,spc,aas[k],spc,dists[j][k]))
end -- if gotit
end -- for j
end -- if 2*htot<tot
print('********** the \'peel-the-onion\' chart above can be used in Excel ('..os.date()..') **********')
print(' ')
--
-- below prepares 2nd chart
if makecsv>1 then
print('Gathering info for the \'AA-stats\' chart below ...')
dstats={} -- will hold dmax davg+drmsd davg davg-drmsd dmin drmsd drms dmax dtopq dmed dbotq dmin for each AA
for j=1,tot do
if j<tot then
dist=dists[j][j+1]
else
dist=dists[j][j-1]
end -- if j
dmin=dist
dmax=dist
dlist={}
dtot=0
for k=1,tot do
if k~=j then
dtot=dtot+1
dist=dists[j][k]
dlist[dtot]=dist
if dist>dmax then
dmax=dist
elseif dist<dmin then
dmin=dist
end -- if dist
end -- if k
end -- for k
if dtot~=tot-1 then
print('ERROR: for j='..j..', dtot='..dtot..'~='..(tot-1)..'=tot-1.')
end -- if dtot
--
-- below gets 5-number summary as in https://en.wikipedia.org/wiki/Five-number_summary
dtopq,dmed,dbotq,davg,drms,drmsd=get5num(dlist,dtot,dmax,dmin,skip5num,sortway)
--
-- gather all the statistics together into dstats
dstats[j]={dmax,davg+drmsd,davg,davg-drmsd,dmin,drmsd,drms,dmax,dtopq,dmed,dbotq,dmin}
end -- for j
--
print(' ')
jnk,js=sortmat(dstats,sortway)
-- now jnk should hold dstats' elements in descending order
-- and js is the sorted order of indices
--
-- below prints 2nd chart with dmax listed in descending order
print('The \'AA-stats\' chart below has one row for each alpha-carbon in the protein.')
print(' Each row lists statistics about the distribution of distances')
print(' from the row\'s alpha-carbon to all other alpha-carbons in the protein:')
print(' dmax=maximum distance, davg=average distance, dmin=minimum distance,')
print(' d=rmsd=rms deviation of all distances from the average distance,')
print(' drms=rms (root-mean-squared) distance (it is used to find d=rmsd),')
print(' dtopq=upper quartile, dmed=median, dbotq=lower quartile distance.')
print(' ')
print('********** the \'AA-stats\' chart below can be used in Excel ('..os.date()..') **********')
tmpstr=string.format('SS%sAA#%sAA%s dmax%sdavg+d%s davg%sdavg-d%s dmin%sd=rmsd%s drms',spc,spc,spc,spc,spc,spc,spc,spc,spc)
print(string.format('%s%s dmax%s dtopq%s dmed%s dbotq%s dmin',tmpstr,spc,spc,spc,spc,spc))
for k=1,tot do
j=js[k]
tmpstr=string.format(' %s%s%3d%s %s',sss[j],spc,j,spc,aas[j])
for m=1,12 do
tmpstr=string.format('%s%s%6.2f',tmpstr,spc,dstats[j][m])
end -- for m
print(tmpstr)
end -- for k
print('********** the \'AA-stats\' chart above can be used in Excel ('..os.date()..') **********')
print(' ')
--
-- below prints 3rd chart
if makecsv>2 then
print('The next chart has one row and one column for each alpha-carbon in the protein.')
print('Each element in the chart is the distance from the alpha-carbon for the element\'s row')
print('to the alpha-carbon for the element\'s column.')
print(' ')
print('********** the \'all-by-all-distances\' chart below can be used in Excel ('..os.date()..') **********')
tmpstr=string.format('SS%sAA#%sAA',spc,spc)
for j=1,tot do
tmpstr=string.format('%s%s%6s',tmpstr,spc,string.format('%s%d%s',sss[j],j,aas[j]))
end -- for j
print(tmpstr)
for j=1,tot do
tmpstr=string.format(' %s%s%3d%s %s',sss[j],spc,j,spc,aas[j])
for k=1,tot do
tmpstr=string.format('%s%s%6.2f',tmpstr,spc,dists[j][k])
end -- for k
print(tmpstr)
end -- for j
print('********** the \'all-by-all-distances\' chart above can be used in Excel ('..os.date()..') **********')
print(' ')
end -- if makecsv>2
end -- if makecsv>1
end -- if makecsv>0
end -- csvcharts()
-- below tests get5num() using examples from:
-- https://en.wikipedia.org/wiki/Quartile (use method 3 results for 1 & 2 below)
-- https://en.wikipedia.org/wiki/Five-number_summary (3)
-- http://www.vias.org/tmdatanaleng/cc_quartile.html (4, but I changed dtopq,dbotq to give method 3 results)
-- http://mathforum.org/library/drmath/view/60969.html (5-9, M&M=method 1, Tukey=method 2)
function test5num(skip5num,sortway)
local dlist,dtot,dmax,dmin,dtopq,dmed,dbotq,davg,drms,drmsd
local j,k,res,tmpstr,resstrs,diff
-- data lists for each example
local dlists={
{6,7,15,36,39,40,41,42,43,47,49},
{7,15,36,39,40,41},
{0,0,1,2,63,61,27,13},
{2,4,7,-20,22,-1,0,-1,7,15,8,4,-4,11,11,12,3,12,18,1},
{1,4,9,16,25,36,49,64,81},
{1,2,3,4,5,6,7,8},
{1,2,3,4,5,6,7,8,9},
{1,2,3,4,5,6,7,8,9,10},
{1,2,3,4,5,6,7,8,9,10,11}}
-- desired values for dmax,dtopq,dmed,dbotq,dmin,drms for each example
-- plus values for dmax,davg+d,davg,davg-d,dmin,d=drmsd for each example
local ress={
{49, 42.75, 40, 20.25, 6, 36.4704,
49, 48.3165, 33.1818, 18.0472, 6, 15.1347},
{41, 40, 37.5, 15, 7, 32.5883,
41, 43.1533, 29.6667, 16.1800, 7, 13.4866},
{63, 44, 7.5, 0.5, 0, 32.7738,
63, 46.1408, 20.8750, -4.3908, 0, 25.2658},
{22, 11.5, 5.5, 0.5, -20, 10.4618,
22, 14.4183, 5.550, -3.3183, -20, 8.8683},
{81, 52.75, 25, 7.75, 1, 41.2755,
81, 58.1410, 31.6667, 5.1924, 1, 26.4743},
{ 8, 6.5, 4.5, 2.5, 1, 5.0498,
8, 6.7913, 4.5, 2.2087, 1, 2.2913},
{ 9, 7.25, 5, 2.75, 1, 5.6273,
9, 7.5820, 5, 2.4180, 1, 2.5820},
{10, 8, 5.5, 3, 1, 6.2048,
10, 8.3723, 5.5, 2.6277, 1, 2.8723},
{11, 8.75, 6, 3.25, 1, 6.7823,
11, 9.1623, 6, 2.8377, 1, 3.1623}}
print('Testing 5-number summary code ('..os.date()..'):')
for j=1,#dlists do
dlist=dlists[j]
dtot=#dlist
dmin=dlist[1]
dmax=dlist[1]
for k=2,#dlist do
if dlist[k]>dmax then
dmax=dlist[k]
elseif dlist[k]<dmin then
dmin=dlist[k]
end -- if dlist
end -- for k
dtopq,dmed,dbotq,davg,drms,drmsd=get5num(dlist,dtot,dmax,dmin,skip5num,sortway)
res={dmax,dtopq,dmed,dbotq,dmin,drms,
dmax,davg+drmsd,davg,davg-drmsd,dmin,drmsd}
resstrs={'dmax',' dtopq','dmed',' dbotq','dmin',' drms',
'dmax','davg+d','davg','davg-d','dmin','d=drmsd'}
for k=1,12 do
if k==1 or k==7 then
tmpstr=''
end -- if k
tmpstr=string.format('%s %s=%6.2f',tmpstr,resstrs[k],res[k])
if k==6 or k==12 then
print(string.format('Example %d gave%s.',j,tmpstr))
end -- if k
end -- for k
for k=1,12 do
diff=res[k]-ress[j][k]
if math.abs(diff)>=0.01 then
print(string.format(' ERROR: %s=%.2f above should be %.2f (a difference of %e).',
resstrs[k],res[k],ress[j][k],diff))
end -- if math.abs(diff)
end -- for k
end -- for j
print('Done testing at '..os.date()..'.')
print(' ')
end -- test5num()
-- below gets 5-number summary as in https://en.wikipedia.org/wiki/Five-number_summary
-- https://en.wikipedia.org/wiki/Median#Easy_explanation_of_the_sample_median
-- https://en.wikipedia.org/wiki/Quartile here Method 3 is used
function get5num(dlist,dtot,dmax,dmin,skip5num,sortway)
local sdlist,hi1,hi2,lo1,lo2,dtopq1,dtopq2,dmed1,dmed2,dbotq1,dbotq2,jnk1,jnk2
local j,dist,davg,drms,drmsd
--
if dtot~=#dlist then
print('ERROR: dtot='..dtot..'~='..(#dlist)..'=#dlist.')
end -- if dtot
--
davg=0
drms=0
for j=1,dtot do
dist=dlist[j]
davg=davg+dist
drms=drms+dist^2
end -- for j
--
-- in the comments below, x is the list of distances in dlist
-- <x> is the average of these distances, or davg
-- N is the # of distances, or dtot
-- <x^2> is the average of the squares of the distances
-- the sums below are over all N values of x
--
davg=davg/dtot
-- now davg = average of dlist's elements
-- average distance: davg=<x>=(sum of x)/N
--
drms=math.sqrt(drms/dtot)
-- now drms = root-mean-squared of dlist's elements
-- rms distance: drms=sqrt[<x^2>]=sqrt[(sum of x^2)/N]
--
drmsd=math.sqrt(drms^2-davg^2)
-- now drmsd = rms deviation of dlist's elements from davg
-- drmsd is like standard deviation but with dtot not dtot-1 in denominator of sqrt
--
-- rms deviation of distances:
-- drmsd=sqrt[(sum of (x-<x>)^2)/N]=root-mean-squared deviation between x and its average, <x>
-- =sqrt[(sum of (x^2-2x<x>+<x>^2))/N]
-- =sqrt[(sum of x^2)/N - 2<x>(sum of x)/N + <x>^2(sum of 1)/N]
-- =sqrt[<x^2> - 2<x><x> + <x>^2 N/N]
-- =sqrt[<x^2> - 2<x>^2 + <x>^2]
-- =sqrt[<x^2>-<x>^2]=sqrt(drms^2-davg^2)
-- standard deviation would use (N-1) in place of N in sqrt[(sum of (x-<x>)^2)/N]
--
if skip5num==2 or (skipit==1 and dtot>=500) then
print('skip5num='..skip5num..' and dtot='..dtot..', so skipping get5num.')
return -1,-1,-1,davg,drms,drmsd -- return crazy values for the quartiles and median of positive #'s
else
-- next line is very rate-limiting if dtot is large like 11476 (puzzle 973 1/17/16)
sdlist,jnk1=sortit(dlist,sortway) -- sdlist is dlist in descending order
if sdlist[1]~=dmax then
print('ERROR: dmax='..dmax..'~='..sdlist[1]..'=sdlist[1].')
end -- if sdlist
if sdlist[dtot]~=dmin then
print('ERROR: dmin='..dmin..'~='..sdlist[dtot]..'=sdlist['..dtot..'].')
end -- if sdlist
hi1,dmed1,lo1=getmedian(sdlist,1)
hi2,dmed2,lo2=getmedian(sdlist,2)
-- methods 1 & 2 above should give dmed1==dmed2
-- but can give different lists for hi1,hi2
-- and different lists for lo1,lo2
jnk1,dtopq1,jnk2=getmedian(hi1,1)
jnk1,dtopq2,jnk2=getmedian(hi2,1)
jnk1,dbotq1,jnk2=getmedian(lo1,1)
jnk1,dbotq2,jnk2=getmedian(lo2,1)
-- method 3 uses the average of methods 1 & 2
return (dtopq1+dtopq2)/2,(dmed1+dmed2)/2,(dbotq1+dbotq2)/2,davg,drms,drmsd
end -- if skip5num
end -- get5num()
-- https://en.wikipedia.org/wiki/Median#Easy_explanation_of_the_sample_median
-- https://en.wikipedia.org/wiki/Quartile here Methods 1 and 2 are used
-- below assumes vals is sorted from hi to lo
function getmedian(vals,method)
local lovals,hivals,med,n,hn,j,midn
n=#vals
hn=n/2
lovals={}
hivals={}
if hn==math.floor(hn) then -- if n is even, hn=n/2=integer
med=(vals[hn]+vals[hn+1])/2
for j=1,hn do
hivals[j]=vals[j]
lovals[j]=vals[j+hn]
-- j+hn rises to 2*hn=2*n/2=n=ok
end -- for j
-- above does not include median in vals, hivals, or lovals
-- and methods 1 & 2 behave the same
else -- n is odd, so midn=(n+1)/2=integer
midn=(n+1)/2
med=vals[midn] -- the median is a datum from the original list in vals
if method==1 then -- do not include median in hivals or lovals (Moore & McCabe, TI-83)
for j=1,midn-1 do
hivals[j]=vals[j]
lovals[j]=vals[j+midn]
-- j+midn rises to (2*midn)-1=(2*(n+1)/2)-1=(n+1)-1=n=ok
end -- for j
elseif method==2 then -- include median in both hivals & lovals (Tukey)
for j=1,midn do
hivals[j]=vals[j]
lovals[j]=vals[j+midn-1]
-- j+midn-1 rises to (2*midn)-1=(2*(n+1)/2)-1=(n+1)-1=n=ok
end -- for j
else
print('ERROR: Method '..method..' is not allowed here.')
end -- if method
end -- if hn
return hivals,med,lovals
end -- getmedian()
-- sortit() w/sortway==1 below was adapted from
-- gary_Minima_Finder_v0.1 code 4/1/14
-- to do like Matlab's sort for a vector x.
-- It should end with x[indx[i]]=s[i] and s[1]>=s[2]>=s[3] etc.
-- so s ends with x's elements in descending order.
--
-- sortway==1 was extremely slow for sorting all 11476 of the
-- distances among 152 alpha-carbons for puzzle 973 on 1/17/16.
-- It looks like it scales slightly below dim^2.
function sortit(x,sortway)
-- sortway==1 does non-table.sort() way
-- sortway==2 does table.sort() way
local s={}
local sx={}
local indx={}
local i,j,k,dim,pflag
dim=#x
pflag=0 -- 0 means no printing so far, 1 means it printed something
if sortway==2 then
-- use table.sort way
for i=1,dim do
sx[i]={}
sx[i][1]=x[i] -- copy x into sx
sx[i][2]=i -- store original indices
end -- for i
table.sort(sx,comp)
-- Above sorts rows of matrix sx in place from max to min
-- using all cols of sx if needed. Having orig index in the
-- last col of sx above makes equal items appear in reverse
-- their original order.
for i=1,dim do
s[i]=sx[i][1] -- retrieve sorted values
indx[i]=sx[i][2] -- retrieve sorted indices
end -- for i
elseif sortway==1 then
-- use non-table.sort way
for i=1,dim do
s[i]=x[i] -- copy x into s
indx[i]=i -- store original indices
end -- for i
k=0
for i=1,dim-1 do
k=k+1
if k==500 then
pflag=1
print(string.format('sortit(%d) is comparing item %5d of %5d with %5d other items (%s).',
sortway,i,dim-1,dim-i,os.date()))
k=0
end -- if k
for j=dim,i+1,-1 do -- i<j always holds here
if s[i]<s[j] then -- use < so get largest element first
s[j],s[i]=s[i],s[j] -- swaps elements i,j
indx[j],indx[i]=indx[i],indx[j]
end -- if s
end -- for j
end -- for i
end -- if sortway
if pflag==1 then
print(' ')
end -- if pflag
return s,indx
end -- sortit()
-- sortmat() below uses the same methods as sortit() above.
-- Its goal is to sort the rows of the matrix x
-- using col 1 1st, then col 2, then col 3, etc, if necessary.
-- It should end with x[indx[i]][j]=s[i][j] and s[1][1]>=s[2][1]>=s[3][1] etc.
-- so s ends with x's rows in descending order.
function sortmat(x,sortway)
-- sortway==1 does non-table.sort() way
-- sortway==2 does table.sort() way
local s={}
local sx={}
local indx={}
local i,j,k,n,nrows,ncols,qflag,swapem,istr,jstr,pflag
nrows=#x
ncols=#(x[nrows])
pflag=0 -- 0 means no printing so far, 1 means it printed something
print(string.format('sortmat(%d) is sorting a matrix with %d rows and %d cols (%s).',sortway,nrows,ncols,os.date()))
pflag=1
if sortway==2 then
-- use table.sort way
for i=1,nrows do
sx[i]={}
for k=1,ncols do
sx[i][k]=x[i][k] -- copy x into sx
end -- for k
sx[i][ncols+1]=i -- store original indices
end -- for i
table.sort(sx,comp)
-- Above sorts rows of matrix sx in place from max to min
-- using all cols of sx if needed. Having orig index in the
-- last col of sx above makes equal items appear in reverse
-- their original order.
for i=1,nrows do
s[i]={}
for k=1,ncols do
s[i][k]=sx[i][k] -- retrieve sorted values
end -- for k
indx[i]=sx[i][ncols+1] -- retrieve sorted indices
end -- for i
elseif sortway==1 then
-- use non-table.sort way
for i=1,nrows do
s[i]={}
for k=1,ncols do
s[i][k]=x[i][k] -- copy x into s
end -- for k
indx[i]=i -- store original indices
end -- for i
n=0
for i=1,nrows-1 do -- i rises from 1
n=n+1
if n==50 then
pflag=1
print(string.format('sortmat(%d) is comparing item %5d of %5d with %5d other items (%s).',
sortway,i,nrows-1,nrows-i,os.date()))
n=0
end -- if n
for j=nrows,i+1,-1 do -- j falls from nrows
-- i<j always holds here
istr=string.format('s[i]=s[%3d]=(',i)
jstr=string.format('s[j]=s[%3d]=(',j)
for k=1,ncols do
istr=string.format('%s %6.2f',istr,s[i][k])
jstr=string.format('%s %6.2f',jstr,s[j][k])
end -- for k
qflag=0
swapem=false -- won't swap i,j below
k=1 -- 1st try to sort using col 1
while qflag==0 do
if s[i][k]<s[j][k] then -- use < so get largest element (here j) first
swapem=true -- will swap i,j below
qflag=1 -- will quit while loop
elseif s[i][k]==s[j][k] then -- if can't sort w/col k, try w/col k+1
k=k+1
if k>ncols then -- have tried all cols, so i,j are same, let em be
qflag=1 -- will quit while loop
end -- if k
else -- s[i][k]>s[j][k] already holds, so no need to swap i,j
qflag=1 -- will quit while loop
end -- if s
end -- while qflag
if k>1 then
pflag=1
print(string.format('sortmat(%d) got to k=%2d of %d sorting s[i]=s[%3d] and s[j]=s[%3d] below (%s):',
sortway,k,ncols,i,j,os.date()))
print(string.format(' %s ).\n %s ).',istr,jstr))
if swapem then
print(' sortmat() will swap the above s[i] and s[j] since s[j] is larger.')
end -- if swapem
end -- if k
if swapem then -- below swaps i,j
indx[j],indx[i]=indx[i],indx[j] -- swaps elements for i,j
for k=1,ncols do
s[j][k],s[i][k]=s[i][k],s[j][k] -- swaps elements for i,j
end -- for k
end -- if swapem
end -- for j
end -- for i
end -- if sortway
if pflag==1 then
print(' ')
end -- if pflag
return s,indx
end -- sortmat()
-- comp() below is used with table.sort() for sortway==2
-- to sort things from highest to lowest.
-- Using comp(si,sj) below put things lowest to highest 1/24/16.
function comp(sj,si)
local qflag,swapem,k,ncols
-- here si=s[i] & sj=s[j], else code is very like in sortmat() for sortway==1
ncols=#si
qflag=0
swapem=false -- won't swap i,j below
k=1 -- 1st try to sort using col 1
while qflag==0 do
if si[k]<sj[k] then -- use < so get largest element (here j) first
swapem=true -- will swap i,j below
qflag=1 -- will quit while loop
elseif si[k]==sj[k] then -- if can't sort w/col k, try w/col k+1
k=k+1
if k>ncols then -- have tried all cols, so i,j are same, let em be
qflag=1 -- will quit while loop
end -- if k
else -- si[k]>sj[k] already holds, so no need to swap i,j
qflag=1 -- will quit while loop
end -- if s
end -- while qflag
return swapem
end -- comp()
-- below eps template is based on 973a19h.eps
-- adapted from eps file made using Matlab
function makeeps1()
print('%!PS-Adobe-2.0 EPSF-1.2')
print('%%Creator: MATLAB, The Mathworks, Inc.')
print('%%BoundingBox: 96 180 501 616')
print('/MathWorks 160 dict begin')
print('/bdef {bind def} bind def')
print('/ldef {load def} bind def')
print('/xstore {exch store} bdef')
print('/c /clip ldef')
print('/cc /concat ldef')
print('/gr /grestore ldef')
print('/gs /gsave ldef')
print('/mt /moveto ldef')
print('/np /newpath ldef')
print('/s {show newpath} bdef')
print('/sg /setgray ldef')
print('/w /setlinewidth ldef')
print('/pgsv () def')
print('/bpage {/pgsv save def} bdef')
print('/epage {pgsv restore} bdef')
print('/bplot /gsave ldef')
print('/eplot {stroke grestore} bdef')
print('/portraitMode 0 def /landscapeMode 1 def /rotateMode 2 def')
print('/dpi2point 0 def')
print('/FontSize 0 def')
print('/FMS {/FontSize xstore findfont [FontSize 0 0 FontSize neg 0 0]')
print(' makefont setfont} bdef')
print('/ISOLatin1Encoding where {pop /WindowsLatin1Encoding 256 array bdef')
print('ISOLatin1Encoding WindowsLatin1Encoding copy pop')
print('/.notdef/.notdef/quotesinglbase/florin/quotedblbase/ellipsis/dagger')
print('/daggerdbl/circumflex/perthousand/Scaron/guilsinglleft/OE/.notdef/.notdef')
print('/.notdef/.notdef/quoteleft/quoteright/quotedblleft/quotedblright/bullet')
print('/endash/emdash/tilde/trademark/scaron/guilsinglright/oe/.notdef/.notdef')
print('/Ydieresis WindowsLatin1Encoding 128 32 getinterval astore pop}')
print('{/WindowsLatin1Encoding StandardEncoding bdef} ifelse')
print('/reencode {exch dup where {pop load} {pop StandardEncoding} ifelse')
print(' exch dup 3 1 roll findfont dup length dict begin')
print(' { 1 index /FID ne {def}{pop pop} ifelse } forall')
print(' /Encoding exch def currentdict end definefont pop} bdef')
print('/isroman {findfont /CharStrings get /Agrave known} bdef')
print('/FMSR {3 1 roll 1 index dup isroman {reencode} {pop pop} ifelse')
print(' exch FMS} bdef')
print('/csm {1 dpi2point div -1 dpi2point div scale neg translate')
print(' dup landscapeMode eq {pop -90 rotate}')
print(' {rotateMode eq {90 rotate} if} ifelse} bdef')
print('/L {lineto stroke} bdef')
print('/MR {4 -2 roll moveto dup 0 exch rlineto exch 0 rlineto')
print(' neg 0 exch rlineto closepath} bdef')
print('/L1i {{currentfile picstr readhexstring pop} image} bdef')
print('currentdict end def')
print('MathWorks begin')
print('bpage')
print('bplot')
print('/dpi2point 12 def')
print('portraitMode 0204 7344 csm')
print('gs 5368 908 133 3406 MR c np')
print('[131 0 0 3404 5369 909] cc')
print('/picstr 1 string def')
end -- makeeps1()
-- below eps template is based on 973a19h.eps
-- adapted from eps file made using Matlab
function makeeps2()
print('5368 908 mt 5368 4313 L')
print('5500 908 mt 5500 4313 L')
print('1874 908 mt 1874 4313 L')
print('1874 908 mt 5279 908 L')
print('1874 4313 mt 5279 4313 L')
print('5279 908 mt 5279 4313 L')
print('eplot')
print('epage')
print('end')
print('showpage')
print('%%EOF')
end -- makeeps2()
function epsmaker(chlist,epsway,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,davg,drmsd,tot,jdim)
local chway,tops,bots,lefts,topys,botys,leftxs,toplo,bothi,leftlo,j10s,j1s,nlines
local j,lnum,tmpstr,tmpstrs,ssstr,aastr,nshades,barys,bary,dcut
if epsway>0 then -- do this routine
print('********** copy lines below to make an eps file ('..os.date()..') **********')
makeeps1() -- set up start of eps file
--
-- below makes middle of eps file
--
if epsway==1 then -- use 5 shades
chway=11
else -- use 10 shades
chway=epsway
end -- if epsway
nshades=#chlist[chway] -- nshades should be 5 or 10
--
-- below shades in the colorbar.
-- oddly 2-digit hex works while 1-digit hex fails here.
-- nevertheless, since hex 0-9A-F give 0-15 in decimal and 255/15=17,
-- the single-hex digits can just be listed twice
-- to get the equivalent 2-digit hex shade.
print(string.format('1 %d 8 [1 0 0 %d 0 0] L1i',nshades,nshades))
-- the 8 in the above line means 8 bits=2^8=256=0-255=00-FF=2-hex digits per shaded box
tmpstr=''
for j=nshades,1,-1 do
tmpstr=string.format('%s%s%s',tmpstr,chlist[chway][j],chlist[chway][j])
end -- for j
print(tmpstr)
print('gr')
--
-- below preps to shade the main part of the plot the 1-digit hex way.
-- would 1-digit hex fail if there were <=10 AA's to be plotted (if jdim<=10) ?
print('gs 1874 908 3406 3406 MR c np')
print('[3404 0 0 3404 1875 909] cc')
print(string.format('/picstr %d string def',jdim))
print(string.format('%d %d 8 [%d 0 0 %d 0 0] L1i',jdim,jdim,jdim,jdim))
-- version 1.1: the 8 in the above line means 8 bits per shaded box
-- (was: the 4 in the above line means 4 bits=2^4=16=0-15=0-F=1-hex digit per shaded box)
--
-- below fills in the graph part with single hexadecimal digits for each value in dds
makechart(chlist,chway,0,0,0,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,jdim,0,true)
print('gr')
--
-- below prepares the text parts for the eps file
tops={} -- up to 7 lines of 60 characters each for above the graph, 1st line is highest
bots={} -- up to 7 lines of 60 characters each for below the graph, 1st line is highest
lefts={} -- up to 7 lines of 48 characters each to left of the graph, 1st line is left-most
tmpstrs={} -- up to 7 lines of 48 characters each, sometimes on top, sometimes on left
tmpstrs[1]=string.format('%s: %s of %s',recipename,user.GetPlayerName(),user.GetGroupName())
tmpstrs[2]=string.format('%s',puzzle.GetName())
tmpstrs[3]=string.format('score=%.3f on %d at %s',trunc3(current.GetScore()),puzzle.GetPuzzleID(),os.date())
tmpstrs[4]=string.format('dist=%.2f (Ca%d-Ca%d) to %.2f (Ca%d-Ca%d)',mindist,minj1,minj2,maxdist,maxj1,maxj2)
tmpstrs[5]=string.format('avg+-rmsdev=%.2f+-%.2f=%.2f to %.2f',davg,drmsd,davg-drmsd,davg+drmsd)
tmpstrs[6]=string.format('epsway=%d w/%d dist ranges for %d of %d AA\'s',epsway,nshades,jdim,tot)
j10s='000000000111111111122222222223333333333444444444455555555556' -- 60 chars wide
j1s='123456789012345678901234567890123456789012345678901234567890' -- 60 chars wide
nlines=math.ceil(tot/60) -- tot=300 will give 5 lines for AAs and 5 lines for SSs
if nlines>5 then -- chop off AAs and SSs beyond 5 lines or with segment #'s > 300
nlines=5
end -- if nlines
if nlines<=2 then -- place all the SS and AA info at the bottom, other info at top
bothi=2*nlines+2
bots[1]=j10s
bots[2]=j1s
toplo=2
for j=1,6 do
tops[toplo+j-1]=tmpstrs[j]
end -- for j
leftlo=8 -- set it above 7 so it won't use the left positions
else -- put the SS info at the top, the AA info at the bottom, and other info to left
bothi=nlines+2
toplo=6-nlines
tops[toplo]=j10s
tops[toplo+1]=j1s
bots[1]=j10s
bots[2]=j1s
leftlo=2
for j=1,6 do
lefts[leftlo+j-1]=tmpstrs[j]
end -- for j
end -- if nlines
for lnum=1,nlines do
j0=(lnum-1)*60
aastr=''
ssstr=''
for j=(j0+1),(j0+60) do
if j<=tot then
aastr=(aastr..aas[j])
ssstr=(ssstr..sss[j])
else
aastr=(aastr..' ')
ssstr=(ssstr..' ')
end -- if j
end -- for j
if nlines<=2 then -- SS and AA info go at bottom
bots[2+lnum]=ssstr
bots[2+lnum+nlines]=aastr
else -- SS info at top, AA info at bottom
tops[toplo+1+lnum]=ssstr
bots[2+lnum]=aastr
end -- if nlines
end -- for lnum
--
-- below are coordinates for different lines of text on the plot
topys={96,216,336,456,576,696,816} -- #'s rise from top to bottom on plot
botys={4476,4596,4716,4836,4956,5076,5196} -- #'s rise from top to bottom
leftxs={1062,1182,1302,1422,1542,1662,1782} -- #'s rise from left to right
barys={4349,4009,3668,3328,2987,2647,2306,1966,1625,1285,944} -- #'s fall bot to top
-- barys are listed in reverse order to correspond better with dcuts & dds
--
-- below adds text to the eps file
print('0 sg')
print('/Courier-Bold /WindowsLatin1Encoding 120 FMSR')
print('6 w')
-- next few lines are to left of plot
for j=leftlo,7 do
print(leftxs[j]..' 4339 mt -90 rotate ('..chopchars(lefts[j],48)..') s 90 rotate')
end -- for j
-- next few lines are above the plot
for j=toplo,7 do
print(string.format('1416 %4d mt (%s) s',topys[j],chopchars(tops[j],60)))
end -- for j
-- next few lines are below the plot
for j=1,bothi do
print('1416 '..botys[j]..' mt ('..chopchars(bots[j],60)..') s')
end -- for j
-- below adds labels and horizontal lines to colorbar
for j=nshades,0,-1 do
if j==0 then dcut='0 '
elseif j==nshades then dcut='Inf'
else
if nshades==5 then -- j=1 to 4, so 2*j=2 to 8
dcut=string.format('%d ',dcuts[2*j])
else -- j=1 to 9, dcuts only has 9 elements
dcut=string.format('%d ',dcuts[j])
end -- if nshades
end -- if j
-- in the above, it is ok to pad the right with extra blanks.
-- we will want just 1st 3 characters in dcut below,
-- and we want the numbers in dcut to be left-justified.
if nshades==5 then -- here j=0 to 5 step 1
bary=barys[2*j+1] -- here 2*j+1=1 to 11 step 2
else -- nshades==10 -- here j=0 to 10 step 1
bary=barys[j+1] -- here j+1=1 to 11 step 1
end -- if nshades
print(string.format('5368 %4d mt 5500 %4d L',bary-36,bary-36))
print(string.format('5520 %4d mt (%3s) s',bary,string.sub(dcut,1,3)))
end -- for j
--
-- above makes middle of eps file
--
makeeps2() -- set up end of eps file
print('********** copy lines above to make an eps file ('..os.date()..') **********')
print(' ')
end -- if epsway
end -- epsmaker()
-- below forces the string str to be dim characters long
function chopchars(str,dim)
local len,res,j,isodd
len=string.len(str) -- len is # of characters in str
if len<dim then -- pad both sides with blank spaces
res=str
isodd=1
for j=1,(dim-len) do
if isodd==0 then -- j is even
res=(' '..res)
else -- j is odd
res=(res..' ')
end -- if j
isodd=1-isodd
end -- for j
elseif len>dim then -- just keep the 1st dim characters of str
res=string.sub(str,1,dim)
else -- just copy str into res
res=str
end -- if len
return res
end -- chopchars()
function makechart(chlist,chway,revflag,spflagi,lbls,sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,jdimi,roflag, doubled )
--
-- version 1.1 -- double up output for eps chart
-- (requires a bit of wrenching)
--
-- add new "doubled" parm, make it optional
--
if doubled == nil then
doubled = false
end
local tmpstr,j100s,j10s,j1s,alleq,ssstr,aastr
local js,j,j1,j2,j100,j10,jrem,jdim,ch,chs,nchs,spflag,spcs,spc
if roflag==1 then -- roflag==1 means use special format for Recipe Output window
jdim=9
spflag=0 -- means no spaces
else
jdim=jdimi
spflag=spflagi -- 0=no spaces, 1=spaces, 2=commas
end -- if roflag
chs={}
nchs = #chlist [ chway ] -- should be 5 or 10
for jj = 1, nchs do
if revflag == 1 then -- do them in reverse order
ch = chlist [ chway ] [ nchs - jj + 1 ]
else -- do them in normal order
ch = chlist [ chway ] [ jj ]
end -- if revflag
-- next part makes sure jj in chs[jj] is 0-9
if nchs == 10 then
chs [ jj - 1 ] = ch
else -- nchs==5
chs [ 2 * jj - 2 ] = ch
chs [ 2 * jj - 1 ] = ch
end -- if nchs
end -- for jj
--
-- lbls==0 all labels off, =1 left labels off, top/bot labels on, =2 all labels on
if lbls==2 then -- left labels are on
j1s=' '
else -- left labels are off
j1s=''
end -- if lbls
if lbls>=1 then -- top/bottom labels are on
j10s=j1s
j100s=j1s
alleq=j1s
ssstr=j1s
aastr=j1s
spcs={'',' ',','}
spc=spcs[spflag+1]
end -- if lbls
js={}
for j=1,jdim do
js[j]=1+round0((j-1)*(tot-1)/(jdim-1)) -- want j=1 to jdim to give js=1 to tot
if lbls>=1 then -- top/bottom labels are on
if j==jdim then
spc=''
end -- if j
j100=trunc0(js[j]/100)
jrem=js[j]-j100*100
j10=trunc0(jrem/10)
j1=jrem-j10*10
j100s=(j100s..j100..spc)
j10s=(j10s..j10..spc)
j1s=(j1s..j1..spc)
alleq=(alleq..'|'..spc)
ssstr=(ssstr..sss[js[j]]..spc)
aastr=(aastr..aas[js[j]]..spc)
end -- if lbls
end -- for j
for j1=jdim,1,-1 do
spcs={'-',' ',','}
spc=spcs[spflag+1]
if roflag==1 then -- do special format for Recipe Output window
tmpstr=string.format('%03d%s',js[j1],spc)
else -- do more general formats
if lbls==2 then -- left labels are on
tmpstr=string.format('%s%03d%s%s',sss[js[j1]],js[j1],aas[js[j1]],spc)
else -- left labels are off
tmpstr=''
end -- if lbls
end -- if roflag
spcs={'',' ',','}
spc=spcs[spflag+1]
for j2=1,jdim do
if j2==jdim then
spc=''
end -- if j2
if doubled then
tmpstr=string.format('%s%s%s%s',tmpstr,chs[dds[js[j1]][js[j2]]],chs[dds[js[j1]][js[j2]]],spc)
else
tmpstr=string.format('%s%s%s',tmpstr,chs[dds[js[j1]][js[j2]]],spc)
end
end -- for j2
if roflag==1 then -- do special format for Recipe Output window
if j1==9 then
tmpstr=string.format('%s %s: for all %d amino acids',tmpstr,recipename,tot)
elseif j1==8 then
tmpstr=string.format('%s min Ca-Ca dist=%6.2f for %3d-%d',tmpstr,mindist,minj1,minj2)
elseif j1==7 then
tmpstr=string.format('%s max Ca-Ca dist=%6.2f for %3d-%d',tmpstr,maxdist,maxj1,maxj2)
elseif j1==6 then
tmpstr=string.format('%s Ca-Ca dist ranges (chway=%d rev=%d):',tmpstr,chway,revflag)
elseif nchs==10 then
if j1==5 then
tmpstr=string.format('%s %s=0 for 0-%-3s %s=1 for %2d-%d',
tmpstr,chs[0],string.format('%d,',dcuts[1]),chs[1],dcuts[1],dcuts[2])
elseif j1==1 then
tmpstr=string.format('%s %s=8 for %2d-%-3s %s=9 for %2d-Inf',
tmpstr,chs[8],dcuts[8],string.format('%d,',dcuts[9]),chs[9],dcuts[9])
else
j=2*(5-j1) -- j1==4 -> j=2, j1==3 -> j=4, j1==2 -> j=6
tmpstr=string.format('%s %s=%d for %2d-%-3s %s=%d for %2d-%d',
tmpstr,chs[j],j,dcuts[j],string.format('%d,',dcuts[j+1]),chs[j+1],j+1,dcuts[j+1],dcuts[j+2])
end -- if j1
else -- nchs==5 should hold
if j1==5 then
tmpstr=string.format('%s %s=0 for 0-%d',tmpstr,chs[0],dcuts[2])
elseif j1==1 then
tmpstr=string.format('%s %s=4 for %2d-Inf',tmpstr,chs[8],dcuts[8])
else
j=(5-j1) -- j1==4 -> j=1, j1==3 -> j=2, j1==2 -> j=3
tmpstr=string.format('%s %s=%d for %2d-%d',tmpstr,chs[2*j],j,dcuts[2*j],dcuts[2*j+2])
end -- if j1
end -- if nchs
end -- if roflag
print ( tmpstr )
end -- for j1
if roflag~=1 then -- do more general formats
if lbls>=1 then -- top/bottom labels are on
print(alleq)
print(ssstr)
print(j100s)
print(j10s)
print(j1s)
print(aastr)
end -- if lbls
end -- if roflag
end -- makechart()
-- below rounds down to 3 decimal places (nearest 0.001)
function trunc3(numi)
local numo=trunc0(numi*1000)/1000
return numo
end -- trunc3()
-- below rounds to 2 decimal places (nearest 0.01)
function round2(numi)
local numo=round0(numi*100)/100
return numo
end -- round2()
-- below rounds to 1 decimal place (nearest 0.1)
function round1(numi)
local numo=round0(numi*10)/10
return numo
end -- round1()
-- below rounds val to the nearest integer to get res
function round0(val)
local res=math.floor(val)
if val-res>=0.5 then
res=res+1
end -- if val
return res
end -- round0()
-- below rounds val down to the nearest integer to get res
function trunc0(val)
local res=math.floor(val)
return res
end -- trunc0()
Main()
