// MACRO CIP
// CIP.spt -- test suite for Jmol CIPChirality.java
// Bob Hanson hansonr@stolaf.edu 4/18/2017 7:44:43 AM
// updated 4/29/2017 1:26:31 PM
// 5/21/2017 4:57:14 PM -- adds AY-236.193,194,195,196 spiro
// 5/27/2017 10:26:49 AM -- finalized - full validation    


VS = [
  id:"VS",
  first: 1,
  last: 300,
  skip: ({})
]



function makeVS() {
  load compounds_3d.sdf
  var n = {*}.model.max
   set echo bottom left
  for (var i = 1; i <= n; i++) {
    frame @i
    echo @i
    refresh
    select visible
    var id = ("00" + i)%-3
    write @{"testSuites/VS/VS_"+id + ".sdf"}
  }

}

function checkDuplicates() {
  load compounds_3d.sdf
  var asmiles = load("slist.txt").lines
  var n = 0;
  for (var s in asmiles) {
    n++;
//    print "" + n + ": " + s
    var a = {*}.find("SMILES/molecule",s)
    if (a.count > 1)
      print "" + n + " " + a.count + "\t" + a[1].model + "\t" + a[2].model + "\t" + s
    else if (a.count == 0)
      print "!!!" + n + " " + a.count  + "\t" + s
    else {
      print "" + n +"\t" + a[1].model + "\t" + s 
    }
  }
}
//checkDuplicates


function labelCIP(n, inverted) {
  if (!n)
    select *
  else
    select atomno=n
  if (inverted)
    invertselected plane xy
  calculate chirality
  labels off
  color labels black
  background labels yellow
  font labels 20
  var a = ("" + n == "false" ? {chirality !=""} : {_C})
  set labelfor @a "%i%[chirality]"
}

function writeMOL(fname, title) {
  select *
  var rs = {chirality != ""}.label("%i%[chirality]").join(" ")
  frame title @{"CIP:" + rs + "; " + title}
  if (!fname.find(".mol"))
    fname += ".mol"
  write coords @fname
  write PNGJ 500 500 @{fname.replace(".mol",".png")}
  molToSdf({}, fname.replace(".mol",".sdf"))
}

sourcedir = "cip" // local version
//sourcedir = "https://www.stolaf.edu/people/hansonr/jmol/cip"

quicksort = true;

// Check a set of Jmol 3D files -- first line includes "CIP:xxxxx;"

makeSDF = false

function loadCIP(info, n) {
  if (0 + info > 0) {
    n = info
    info = thisSDF
  }
  thisSDF = info
  zap
  if (info.sdffile) {
      load @{info.sdffile} @i
      print _M.molData
  }
}

function checkCIP(info, n) {
  if (0 + info > 0) {
    n = info
    info = thisSDF
  }
  thisSDF = info
  var first = (n ? n : info.first)
  var last = (n ? n : info.last)
  var bb
  for (var i = first; i <= last; i++) {
    if (info.skip & i) {
      print "SKIPPING MODEL " + info.id + " " + i
      continue
    }
    if (makeSDF) {
      if (info.sdffile) {
        load @{info.sdffile} @i
        print _M.molData
        bb = (_M.molData ? _M.molData.ref_pname_pin : "")
      } else {
        bb = ""
        if (_modelTitle.find("BB:"))
          bb = _modelTitle.split("BB:")[2]
      }
    }
    var id = info.id + "_" + (("000" + i)%-3)
    var f = sourceDir + "/testSuites/"+info.id+"/" + id + ".sdf";
    if (checkrs(f, 1, "*") && makeSDF) {
      var molData = _M.molData
      var note = (molData ? molData.jmol_note : "")
      molData = {}
      if (bb)
        molData.bb_ref = bb
      molData.jmol_id = id
      molToSdf(molData)
    }
  }
}

function getJmolCIP() {
  return {chirality != ""}.label("%i%[chirality]").join(" ")
}

function molToSdf(molData, fname) {
  if (!molData) molData = {}
  select *
  if (!fname)
    fname = _modelFile.replace(".mol",".sdf");
  calculate chirality
  var rs = getJmolCIP()
  molData.jmol_cip = rs
  molData.jmol_smiles = {*}.find("SMILES/noaromatic")
  molData.nci_smiles = show("chemical smiles")
  // for first line of mol/sdf file
  var title = (molData.jmol_id ? 'ID:' + molData.jmol_id + ";" : "") 
    + 'CIP:' + rs + ";"
    + (molData.bb_ref ? bb ? 'BB:' + molData.bb_ref + ";" : "")    
  frame title @title 
  model property "molData" molData
  write coord @fname
  print molData
  print ">>\t" + molData.jmol_id + "\t" + rs + "\t" + molData.nci_smiles
}

// check a set of files in a directory based on a file listing with annotations

function checkRdir(name, type) {
  var x = load(sourceDir + "/" + name + ".txt").lines
  for (var f in x) {
    f = f.trim();
    if (f == "#QUIT") break
    if (f == "#stop") exit
    if (!f || f.find("#") == 1) continue
    if (f.find("$") != 1) 
      f = sourceDir + "/" + name + "/" + f      
    if (type)
      checkRS(f, 1, type)
    else
      checkRS(f)
  }
}
  
function showRS(prettyN, justChiralLabels) {
  
  select !_H
  color labels black
  background labels yellow
  set labelfor {_C} "%a"
  refresh
  if (prettyN)
    select atomno=prettyN
  else
    select *
  calculate chirality; 
  label ""
  // leave s(5) [and r(5)] because it shows the difference between S and s
  set labelfor {chirality != ''} "%a %[chirality]"
  set labelfor {chirality != '' && ciprule != '1a'} "%a %[chirality](%[ciprule])"
  select *
  font label 20
  set labelfront
  set labeloffset 10 10
  
  if (prettyN) {

    display remove _H

    select *
    set labelfront
    set labeloffset 0 0
    background labels none
    if (justChiralLabels) {
      set labelfor {chirality == ''} ""
      select {chirality == ''}
      spacefill off
    }

    select atomno=prettyN
    font label 35
    spacefill 0.3
    set labelfront
    display add connected(atomno=prettyN)
  }
}

nOK = 0;

thisN = 0;
stopped = false;

function checkRS(fname, n, key, isQuiet) {
  if (_argCount == 2 && 0 + n != n) {
    key = n;
    n = 0;
  }
  if (!n)
    n = 1;
  var doCheck = !!key
  fname = fname.replace("://",":/")
  if (fname.find("//")) {
    fname = fname.split("//")
    key = fname[2]
    fname = fname[1]
    doCheck = key
  }
  fname = fname.replace(":/", "://")
  if (thisN && (nOK + 1) != thisN) {
   nOK++;
   return false;
  }
  print "loading " + (nOK + 1) + ": " + fname + " " + n + " k=" +  key + " c=" + docheck + " q=" + quicksort
  set useMinimizationThread false
  load @fname @n filter "2D"
  if (!{_H})  calculate hydrogens
  
  
  if (docheck && quicksort) {
    var t = now();
    calculate chirality
    totalms += now(t);
    select (chirality != "")
  } else {  
    set labelfor {_C} "%[atomname]"
    rotate best
    color labels black
    background label yellow
    refresh
    calculate chirality
    select (chirality != "")
    set labelfor {_C} ""  
    set labelfor {selected} "%[atomname] %[chirality]"
    set labelfor {selected && ciprule != '1a'} "%[atomname] %[chirality](%[ciprule])"
    set labelfront
    set labeloffset 10 10
  }
  
  var rs = {selected}.label("%[chirality]").join("")
  
  var isok = true
  var isjmol = false
  if (doCheck) {
    var ref = ""
    var refvs = ""
    if (docheck != key) {
      if (_M.molData) {
        ref = _M.molData["jmol_cip"];
      }
      if (ref) {
        isjmol = true 
      } else if (_modelTitle.find("CIP:")) {
        ref = _modelTitle.split("CIP:")[2].split(";")[1]
        isjmol = true
      } else if (_M.molData) {
        // accepts data in various forms
        ref = _M.molData.cip_labels
        if (ref) {
          isVS = true;
        } else {
          ref = _M.molData["chiral_atoms"].replace("\n","").replace(" ","") + _M.molData["stereounits"] + _M.molData["stereo"]
        }      
      }
    }
    if (ref) {
      rs = ""
      if (isjmol || isVS){
        rs = {selected}.label("%i%[chirality]").join(" ")
      } else if (ref.find("focus:")) {
        // MY64
        var tokens = ref.split("focus:");
        ref = ""
        var refs = []
        for (var focus in tokens) {
          if (!focus) continue
          var items = focus.split("descriptor:")
          var atoms = items[1].split(true)
          rs = items[2].trim()
          var na = atoms.length
          refs[1 + atoms[1]] = rs // first
          refs[1 + atoms[0]] = rs // last
        }
        var na = refs.length
        for (var j = 1; j <= na; j++) {
          if (refs[j])
            ref += "" + j + refs[j]
        }
        rs = ""
      } else if (ref.find(",")) {
        rs = {!_H}.chirality.join(",").trim(",")
        ref = ref.trim(",")
      } else if (_M.molData["chiral_atoms"]) {
        rs = ""
      } else {
        rs = {!_H}.chirality.join("")
      }
      // switch to BB2013 preferred name format 
      key = ref.replace("Sa","P").replace("Ra","M").replace("sa","p").replace("ra","m");
      if (!rs)
        rs = {selected}.label("%i%[chirality]").join("")
    } else if (key == "*") {
      key = ""
    }
    
    if (key.find(",") && !rs.find(","))
      key = key.replace(",","")
      
    // oddly, "like" in Jmol means "isExactly"; just is
    if (key like rs) {
      nOK++
      print _M.molData
      print "" + nOK + " OK\t" + fname + "\t" + rs
    } else {
      refresh 
      var ans = "yes"
      isOK = false;
      if (isQuiet) {
        print ',"' + n + '":"' + rs + '" // ' + key  
      } else {
        var s = "" + n + "??\t" + fname + "\t // found " + rs + "; should be " + key
        print s
//        ans = "yes"
        ans = prompt(s.replace("\t"," ") + " \n\n continue?", "yes")
      }
      if (ans != "yes") quit 
    }
  } else {
      print fname + "\t" + rs
  }
  select *
  refresh
  return isOK
}

///// older code for extracting from SDF files

function findSDF(info, what) {
  if (!what) {
    what = info
    info = ""
  }
  if (!info) info = thisSDF
  load @{info.sdffile}
  var n = {*}.model.max
  var data = getProperty("auxiliaryInfo")
  var found = []
  var ivals = []
  for (var i = 1; i <= n; i++) {
    var a = data.models[i].molData
    if (a.find(what)) {
      print "" + i + ": " + a
      found.push(a)
      ivals.push(i)
    }
  }
  print "" + found.length + " found for " + what + ": " + format("json", ivals) 
}

function bhtest(info, n) {
  // bhtest 3
  // bhtest MV64 3
  if (info) {
    if (!n) {
      n = info
      info = ""
    }
  }
  checkCIP((info ? info : thisSDF),n)
  label %a;background labels yellow
  refresh
  select *
  //set loglevel 6
  print {*}.chirality.join(",").trim(",")
  //print {*}.chirality.pivot
  print {*}.find("SMILES")
  showRS
}



/////////////////////////////////////////////////////////////////////
//
// structures from a multi-model SDF file that may have stereochemistry annotations
//
// @3D means there is a file testSuites/3D/test<n>.mol
//
// Test Suite AY-236 (Andrey Yerin, ACD/Labs, Moscow, Russia)
//
// BB errata page: http://www.chem.qmul.ac.uk/iupac/bibliog/BBerrors.html

AY236 = [
  id:"AY236",
  sdffile:"c:/temp/P-9-examples-stereo.sdf",
  first: 1,
  last:236,
  skip:({27 229}) || // ignore -- BB has E/Z only; missing chirality
         ({95 96 98 99 100 101 102 103 104 108 109 110 111 112 200}) || // trigonal planar, square planar, or hypervalent
         ({212 213})|| // chiral conformation 1,4-benzene in a ring
         ({38}) // ignoring -- 38: Mancude for cyclopentadienyl anion
]

//
// Test Suite MV-64 (Mikko Vainio,  �bo Akademi University, Turku, Finland)
//

MV64 = [
  id:"MV64",
  sdffile:"testSuites/MV64/stereo_test_cases.sdf",
  first: 1,
  last: 64,
  skip: ({34 35 36 37}) // pentacoordinate P 
]

//
// Test Suite MV-116 (Mikko Vainio,  �bo Akademi University, Turku, Finland)
//
// ZIP file cleaned of redundant 2D files, hypervalent models, and undefined stereochemistry, planar chirality

MV116 = [
  id:"MV116",
  zipfile:"testSuites/MV116/cip_examples.zip",
  first: 1,
  last: 80,
  skip: ({})
]

BH64 = [
  id:"BH64",
  first: 1,
  last: 83,
  skip: ({})
]

JM = [
  id:"JM",
  first: 1,
  last: 9,
  skip: ({})
]

L = [
  id:"L",
  first: 35,
  last: 87,
  skip: ({36:55 57:84 86})
]




thisSDF = AY236

print "checkCIP loaded. presumes local AY236 in the cip subdirectoryy. Examples:"
print "\tbhtest 53"
print "\tbhtest AY236 132"
print "\twriteMol xxx.mol"
print "\tshowRS"
print "\tlabelCIP"
print "\tlabelCIP false"

print "\tcheckCIP(AY236)"
print "\tcheckCIP(MV116)"
print "\tcheckCIP(MV64)"
print "\tcheckCIP(BH64)"
print "\tcheckCIP(VS)"




/////////////////////////////////////////////////////////////////////

var ans = prompt("Do you want to run all tests?", "no")
if (ans == "null" || ans == "no") quit

// full test run

totalms = 0

set testflag1 true // skips data tracking
set ciprule6full false // would enable rr

if (ans == "VS") {
  checkCIP(VS)
} else {

  checkCIP(BH64)
  checkCIP(AY236)
  checkCIP(MV116)
  checkCIP(MV64)
  checkCIP(JM)
  checkCIP(L)
  
}  
set testflag2 false

//checkRdir("bytype/more", "*")
//checkRdir("bytype/centres", "*")

print "DONE: " + nOK + " ms = " + totalms 




/**

// individual checks 

// checkRS("$glucose", "RSRR")
checkRS("$(2S,3R)-2,3-oxiranediol", "SR")
checkRS("$(S)-2-butanol", "S")
checkRS("$(R)-2-butanol", "R")
checkRS("$(2S,3R)-2,3-butanediol", "SR")
checkRS("$(2S,3S)-2,3-butanediol", "SS")
checkRS("$(2R,3R)-2,3-butanediol", "RR")
checkRS("$(2R,3S)-2,3-butanediol", "RS")
checkRS("$1,4-dimethylcyclohexane", "")
checkRS("$cholesterol", "RRSSSRSR") // (3S,8S,9S,10R,13R,14S,17R) and sidechain (R)
checkRS("==ta1", "SSRSRSSRSRS") // taxol (1S,2S,3R,4S,7R,9S,10S,12R,15S) and sidechain (2R,3S)
checkRS("testSuites/VS/VS_209.mol")


load $OC(/[16OH])=C(/O)[18OH]    
showrs
select _F or _O
background labels blue
color labels white
label %[mass]