context: using System;
  using System.Collections.Generic;
  using System.Linq;
  using System.Text;
  using uQlustCore;
  using uQlustCore.PDB;
  
  
  namespace uQlustCore.Distance
  {
      public class MaxSub:Rmsd
      {
  
          public MaxSub(DCDFile dcd, string alignFile, bool flag, string refJuryProfile = null)
              : base(dcd, alignFile, flag, PDBMODE.ONLY_CA,refJuryProfile)
          {
              order = true;
              maxSimilarity = 0.0;
          }
          public MaxSub(string dirName, string alignFile, bool flag, string refJuryProfile = null)
              : base(dirName, alignFile, flag, PDBMODE.ONLY_CA,refJuryProfile)
          {
              this.dirName = dirName;
              this.alignFile = alignFile;
              this.flag = flag;
              this.refJuryProfile = refJuryProfile;
  
              order = true;
              maxSimilarity = 0.0;
          }
          public MaxSub(List<string> fileNames, string alignFile, bool flag, string refJuryProfile = null)
              : base(fileNames, alignFile, flag, PDBMODE.ONLY_CA, refJuryProfile)
          {
              
              order = true;
              maxSimilarity = 0.0;
          }
          public override string ToString()
          {
              return "MaxSub";
          }
          public override bool SimilarityThreshold(float threshold, float dist)
          {
              if (dist < threshold)
                  return true;
              return false;
  
          }
          public override double GetRealValue(double v)
          {
              return (100 - v) / 100;
          }
          protected KeyValuePair<List<int>,float[,]> FindLongestSegment(int seedL,float threshold, float [,] mol1,float [,] mol2)
          {
              List<int> bestList = null;
              List<int> pList = new List<int>();
              List<int> cList = null;
              float[,] remTransMol = null;
              float[,] copyMol2 = null;
              copyMol2 = new float[mol2.GetLength(0), mol2.GetLength(1)];
              remTransMol = new float[mol1.GetLength(0), mol1.GetLength(1)];
              for (int i = 0; i < mol1.GetLength(0) - seedL; i++)
              {
                  pList.Clear();
                  for (int j = 0; j < seedL; j++)
                      pList.Add(i + j);
  
                  KeyValuePair<List<int>, float[,]> extRes = Extend(pList, bestList, mol1, mol2, copyMol2,threshold);
                  cList = extRes.Key;
  
                  if (cList == null)
                      continue;
  
                  if (bestList==null || cList.Count > bestList.Count)
                  {
                      bestList = new List<int>(cList);
                      Buffer.BlockCopy(extRes.Value, 0, remTransMol, 0, copyMol2.Length * sizeof(float));
                  }
                  if (cList.Count == mol1.GetLength(0))
                      break;
  
  
              }
              KeyValuePair<List<int>, float[,]> res = new KeyValuePair<List<int>, float[,]>(bestList, remTransMol);
  
              return res;
          }
  
          public override int GetDistance(string refStructure, string modelStructure)
          {
              int seedL = 4;
              float maxDist = 3.5f;
  
              if (!pdbs.molDic.ContainsKey(refStructure) || !pdbs.molDic.ContainsKey(modelStructure))
                  return errorValue;
  
              
              posMOL locPosMol = Optimization.PrepareData(pdbs.molDic[refStructure], pdbs.molDic[modelStructure]);
  
              KeyValuePair<List<int>,float[,]> seg = FindLongestSegment(seedL, maxDist, locPosMol.posmol1, locPosMol.posmol2);
              if ( seg.Key == null)
                  return 0;
  
              float sum = 0;
              float mDist = maxDist * maxDist;
              float rmsd = 0;
                  foreach(var item in seg.Key)
                  {
                      float dist = 0;
                      dist = calcDist(locPosMol.posmol1, seg.Value, item);
                      rmsd += dist;
                      dist *= dist;
                      sum += 100 / (1 + dist / mDist);
                  }
              rmsd /= seg.Key.Count;
              return 100 - (int)((sum / pdbs.molDic[refStructure].molLength));//locPosMol.posmol1.GetLength(0)));
          }
  
          //Warning MaxSub has been reversed 0 is the best one, 1 is the worst
          /*public override int GetDistance(string refStructure, string modelStructure)
          {
              int seedL = 4;
              int maxS = 0;
              float maxDist = 3.5f;
              //float maxDist = 1.5f;
              List<int> pList = new List<int>();
              List<int> cList = null;
              float[,] remMol = null;
              float[,] remTransMol = null;
              float[,] copyMol2 = null;
              List<int> bestList = null;
  
           //   TestRmsd(refStructure, modelStructure);
              if (!pdbs.molDic.ContainsKey(refStructure) || !pdbs.molDic.ContainsKey(modelStructure))
                  return int.MaxValue;
  
              posMOL locPosMol=Optimization.PrepareData(pdbs.molDic[refStructure], pdbs.molDic[modelStructure]);
              copyMol2 = new float[locPosMol.posmol2.GetLength(0), locPosMol.posmol2.GetLength(1)];
  
              remMol = new float[locPosMol.posmol1.GetLength(0), locPosMol.posmol1.GetLength(1)];
              remTransMol = new float[locPosMol.posmol1.GetLength(0), locPosMol.posmol1.GetLength(1)];
  
              //opt.CenterMols(opt.posMol1, opt.posMol2);
              if (bestList != null)
                  bestList.Clear();
              for (int i = 0; i < locPosMol.posmol1.GetLength(0) - seedL; i++)
              {
                  pList.Clear();
                  for (int j = 0; j < seedL; j++)
                      pList.Add(i + j);
  
                  KeyValuePair<List<int>, float[,]> extRes = Extend(pList, bestList,locPosMol.posmol1, locPosMol.posmol2, copyMol2, maxDist);
                  cList = extRes.Key;
  
                  if (cList == null)
                      continue;
  
                  if (cList.Count > maxS)
                  {
                      bestList = new List<int>(cList);
                      maxS = cList.Count;
               //       Buffer.BlockCopy(copyMol1, 0, remMol, 0, copyMol2.Length * sizeof(float));
                      Buffer.BlockCopy(extRes.Value, 0, remTransMol, 0, copyMol2.Length * sizeof(float));                    
                  }
                  if (cList.Count == locPosMol.posmol1.GetLength(0))
                      break;
  
  
              }
              if (remTransMol == null || bestList==null || bestList.Count==0)
                  return int.MaxValue;
  			
  			
  			
              float sum = 0;
              float mDist = maxDist * maxDist;
  			float rmsd=0;
              foreach (var item in bestList)
              {
                  float dist = 0;
                  dist = calcDist(locPosMol.posmol1, remTransMol, item);
  				rmsd+=dist;
                  dist *= dist;
                  sum += 100/(1 + dist/ mDist);
              }
  			rmsd/=bestList.Count;
              return 100 - (int)((sum / locPosMol.posmol1.GetLength(0)));
          }*/
          private float [,] Rotate(List<int> indexList, float[,] pMol1, float[,] pMol2,float [,]copyMol)
          {
              float[,] indexMol1;
              float[,] indexMol2;
              float[,] transMatrix = null;
              float[,] rotMol = null;
  			
  			
  
  			if (indexList.Count <= 2)
                  return null;
  
  			
              indexMol1 = new float[indexList.Count, 3];
              indexMol2 = new float[indexList.Count, 3];
              for (int j = 0; j < indexList.Count; j++)
              {
                  for (int l = 0; l < 3; l++)
                  {
                      indexMol1[j, l] = pMol1[indexList[j], l];
                      indexMol2[j, l] = pMol2[indexList[j], l];
                  }
              }
              float[] center1 = new float[3];
              float[] center2 = new float[3];
              //opt.CenterMols(indexMol1, indexMol2);
              Optimization.CenterMol(indexMol1, center1);
              Optimization.TransVec(indexMol1, center1);
              Optimization.CenterMol(indexMol2, center2);
              Optimization.TransVec(indexMol2, center2);
  
  
gt:             transMatrix = Optimization.TransMatrix(indexMol1,indexMol2);