sgss.c

/* $Id: sgss.c 851 2002-02-09 05:39:17Z wdelano $ */

/* The source code contained in this file is            */
/* Copyright (C) 1994-2000 by Ralf W. Grosse-Kunstleve. */
/* Please see the LICENSE file for more information.    */

#ifdef _PYMOL_WIN32
#include"os_predef.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


#undef SG_GLOBAL
#include "sglite.h"
#include "sgconst.h"


#define  mDiscrGr  8


typedef struct
  {
    int  P[3];
    int  Z[3];
  }
  T_DiscrGr;


int Is_ss(const T_ssVM *ssVM, int n_ssVM, int h, int k, int l)
{
  int  i_ssVM, u;

  range1(i_ssVM, n_ssVM)
  {
    u =  ssVM[i_ssVM].V[0] * h;
    u += ssVM[i_ssVM].V[1] * k;
    u += ssVM[i_ssVM].V[2] * l;

    if (ssVM[i_ssVM].M) {
      if (u % ssVM[i_ssVM].M) return 0; }
    else {
      if (u)                  return 0; }
  }

  return 1;
}


void Set_uvw(const T_ssVM *ssVM, int n_ssVM, int h, int k, int l, int uvw[3])
{
  int  i_ssVM, u;

  range1(i_ssVM, n_ssVM)
  {
    u =  ssVM[i_ssVM].V[0] * h;
    u += ssVM[i_ssVM].V[1] * k;
    u += ssVM[i_ssVM].V[2] * l;

    if (ssVM[i_ssVM].M) u %= ssVM[i_ssVM].M;

    uvw[i_ssVM] = u;
  }
}


/* sgss.c:SetAnyIxGen() and sginfo2/metric.c:SetIxGen() should be consolidated.
 */
static int SetAnyIxGen(const T_SgOps *SgOps, int PG_MGC, int IxGen[2])
{
  int     nGen, iLSMx, i;
  int     PrincipalProperOrder;
  T_RMxI  PrincipalRI[1], RI[1];

  rangei(2) IxGen[i] = -1;

  if (PG_MGC == MGC_Undefined) PG_MGC = GetPG(SgOps);
  if (PG_MGC == MGC_Unknown) return IE(-1);

  nGen = 0;
  PrincipalProperOrder = 0;

  switch (ixXS(PG_MGC))
  {
    case XS_Triclinic:
      if (SgOps->fInv == 1) {
        IxGen[0] = 0;
        nGen = 1;
      }
      break;

    case XS_Monoclinic:
      IxGen[0] = 1;
      nGen = 1;
      break;

    case XS_Orthorhombic:
      IxGen[0] = 1;
      IxGen[1] = 2;
      nGen = 2;
      break;

    case XS_Tetragonal:
                                  PrincipalProperOrder = 4;
    case XS_Trigonal:
      if (! PrincipalProperOrder) PrincipalProperOrder = 3;
    case XS_Hexagonal:
      if (! PrincipalProperOrder) PrincipalProperOrder = 6;

      range2(iLSMx, 1, SgOps->nSMx)
      {
        if (SetRotMxInfo(SgOps->SMx[iLSMx].s.R, PrincipalRI) == 0)
          return IE(-1);

        if (abs(PrincipalRI->Rtype) == PrincipalProperOrder) {
          if (PrincipalRI->SenseOfRotation > 0) {
            IxGen[0] = iLSMx;
            nGen++;
            break;
          }
        }
      }

      if (nGen < 1) return IE(-1);

      range2(iLSMx, 1, SgOps->nSMx)
      {
        if (iLSMx == IxGen[0]) continue;

        if (SetRotMxInfo(SgOps->SMx[iLSMx].s.R, RI) == 0)
          return IE(-1);

        if (   abs(RI->Rtype) == 2
            && MemCmp(RI->EV, PrincipalRI->EV, 3) != 0) {
          IxGen[1] = iLSMx;
          nGen++;
          break;
        }
      }

      break;

    case XS_Cubic:
      range2(iLSMx, 1, SgOps->nSMx)
      {
        if (SetRotMxInfo(SgOps->SMx[iLSMx].s.R, RI) == 0)
          return IE(-1);

        if      (   abs(RI->Rtype) == 3
                 && RI->SenseOfRotation > 0
                 && IxGen[0] < 0) {
                    IxGen[0] = iLSMx;
          nGen++;
        }
        else if (   abs(RI->Rtype) == SgOps->nSMx / 6
                 && RI->SenseOfRotation >= 0
                 && IxGen[1] < 0) {
                    IxGen[1] = iLSMx;
          nGen++;
        }
      }

      if (nGen != 2) return IE(-1);

      break;

    default:
      return IE(-1);
  }

  return nGen;
}


static int CmpOLen2(const int a[3], const int b[3])
{
  int OLen2a, OLen2b, i;

  OLen2a = 0; rangei(3) OLen2a += a[i] * a[i];
  OLen2b = 0; rangei(3) OLen2b += b[i] * b[i];

  if (OLen2a < OLen2b) return -1;
  if (OLen2a > OLen2b) return  1;

  return CmpiVect(a, b, 3);
}


static int ConstructGenRmI(const T_SgOps *SgOps, const T_RTMx Z2PCBMx[2],
                           const int IxGen[2], int nGen,
                           int GenRmI[3 * 3 * 3])
{
  int  nrGenRmI;
  int  iRmI, iGen, i;

      nrGenRmI = (SgOps->fInv - 1 + nGen) * 3;
  if (nrGenRmI > 9) return IE(-1);

  iRmI = 0;

  if (SgOps->fInv == 2) {
    SetRminusI(SgOps->SMx[0].s.R, GenRmI, 1);
    iRmI++;
  }

  if (Z2PCBMx == NULL) {
    range1(iGen, nGen) {
      SetRminusI(SgOps->SMx[IxGen[iGen]].s.R, &GenRmI[iRmI * 9], 0);
      iRmI++;
    }
  }
  else {
    range1(iGen, nGen) {
      if (CB_RMx(&GenRmI[iRmI * 9],
                 Z2PCBMx[0].s.R,
                 SgOps->SMx[IxGen[iGen]].s.R,
                 Z2PCBMx[1].s.R) != 0) return -1;
      range3(i, 0, 9, 4) GenRmI[iRmI * 9 + i] -= 1;
      iRmI++;
    }
  }

  if (iRmI * 3 != nrGenRmI) return IE(-1);

  return nrGenRmI;
}


static int GetContNullSpace(const T_SgOps *SgOps, int IxGen[2], int nGen,
                            T_ssVM ssVM[3])
{
  int  GenRmI[3 * 3 * 3];
  int  nrGenRmI, RankGenRmI;
  int  nIndep, iIndep, IxIndep[3], Sol[4][3];
  int  n_ssVM;

      nrGenRmI = ConstructGenRmI(SgOps, NULL, IxGen, nGen, GenRmI);
  if (nrGenRmI < 0) return IE(-1);

      RankGenRmI = iRowEchelonFormT(GenRmI, nrGenRmI, 3, NULL, 0);
  if (RankGenRmI < 0 || RankGenRmI > 3)
    return IE(-1);

  n_ssVM = 3 - RankGenRmI;

      nIndep = iRESetIxIndep(GenRmI, RankGenRmI, 3, IxIndep, 3);
  if (nIndep < 0) return IE(-1);

  if (nIndep != 2)
  {
    range1(iIndep, nIndep)
    {
      ssVM[iIndep].V[IxIndep[iIndep]] = 1;

      if (iREBacksubst(GenRmI, NULL, RankGenRmI, 3,
                       ssVM[iIndep].V, NULL) < 1)
        return IE(-1);

      ssVM[iIndep].M = 0;
    }
  }
  else
  {
    if (SolveHomRE1(GenRmI, IxIndep, Sol) != 0) return -1;

    qsort((void *) Sol, 4, sizeof (*Sol),
          (int (*)(const void *, const void *)) CmpOLen2);

    range1(iIndep, 2) {
      MemCpy(ssVM[iIndep].V, Sol[iIndep], 3);
      ssVM[iIndep].M = 0;
    }
  }

  return n_ssVM;
}


static int nDLoopStep(int *i, int n, int Low, int High)
{
  int  p, l;

  p = l = n - 1;

  for (; p >= 0;)
  {
        i[p]++;
    if (i[p] > High)
      p--;
    else if (p < l)
      i[++p] = Low - 1;
    else
      return 1;
  }

  return 0;
}


static void UpdateBestZ(int OrigZf[mDiscrGr][3], int nDiscrGr,
                        int BestZf[mDiscrGr][3], int BestM[mDiscrGr],
                        int BestZc[mDiscrGr][3],
                        int Shift[3], int LTBF)
{
  int  iDG, c, i;
  int  Zf[3], M;
  int  Zc[3];

  range2(iDG, 1, nDiscrGr)
  {
    rangei(3) Zf[i] = iModPositive(Shift[i] + OrigZf[iDG][i], LTBF);
    MemCpy(Zc, Zf, 3);
    M = CancelBFGCD(Zc, 3, LTBF);

    rangei(3) {
      if (Zf[i]) {
        c = CmpOLen2(BestZc[iDG], Zc);
        if (c > 0 || (c == 0 && BestM[iDG] > M)) {
          MemCpy(BestZf[iDG], Zf, 3);
          MemCpy(BestZc[iDG], Zc, 3);
          BestM[iDG] = M;
        }
        break;
      }
    }
  }
}


static int BestVect(const T_SgOps *SgOps,
                    const T_ssVM *ssVM, int n_ssVM,
                    int DTBF, T_DiscrGr DiscrGr[mDiscrGr], int nDiscrGr)
{
  int  iLTr, i_ssVM, iDG, gcd, i, j;
  int  fGrd, LTBF, f[2];
  int  OrigZf[mDiscrGr][3];
  int  BestZf[mDiscrGr][3], BestM[mDiscrGr];
  int  BestZc[mDiscrGr][3];
  int  LTr[3][3];

  fGrd = 1;
  LTBF = 1;

  range2(iDG, 1, nDiscrGr) {
    rangei(3) {
      gcd = iGCD(DiscrGr[iDG].Z[i], DTBF * CRBF);
      LTBF = iLCM(LTBF, (DTBF * CRBF) / gcd);
    }
  }

  range2(iLTr, 1, SgOps->nLTr) {
    rangei(3) {
      gcd = iGCD(SgOps->LTr[iLTr].v[i], STBF);
      LTBF = iLCM(LTBF, STBF / gcd);
    }
  }

  range1(i_ssVM, n_ssVM)
    rangei(3) fGrd = iLCM(fGrd, ssVM[i_ssVM].V[i]);

  LTBF *= fGrd;

  if (LTBF > 6) LTBF = iLCM(LTBF,  6);
  else          LTBF = iLCM(LTBF, 12);

  if (SgOps->nLTr == 1 && n_ssVM == 0) return 0;

  range2(iDG, 1, nDiscrGr) {
    if (ChangeBaseFactor(DiscrGr[iDG].Z, DTBF * CRBF,
                            OrigZf[iDG],        LTBF, 3) != 0)
      return IE(-1);
    rangei(3) OrigZf[iDG][i] = iModPositive(OrigZf[iDG][i], LTBF);
    MemCpy(BestZf[iDG], OrigZf[iDG], 3);
    MemCpy(BestZc[iDG], OrigZf[iDG], 3);
    BestM[iDG] = CancelBFGCD(BestZc[iDG], 3, LTBF);
  }

  if (n_ssVM > 2) return IE(-1);

  range1(iLTr, SgOps->nLTr)
  {
    if (ChangeBaseFactor(SgOps->LTr[iLTr].v, STBF, LTr[0], LTBF, 3) != 0)
      return IE(-1);

    rangei(n_ssVM) f[i] = 0;

    do
    {
      rangei(n_ssVM)
        range1(j, 3) LTr[i + 1][j] = LTr[i][j] + f[i] * ssVM[i].V[j];

      UpdateBestZ(OrigZf, nDiscrGr, BestZf, BestM, BestZc,
                  LTr[n_ssVM], LTBF);
    }
    while (nDLoopStep(f, n_ssVM, 0, LTBF - 1));
  }

  range2(iDG, 1, nDiscrGr)
    if (ChangeBaseFactor(BestZf[iDG],           LTBF,
                         DiscrGr[iDG].Z, DTBF * CRBF, 3) != 0)
      return IE(-1);

  return 0;
}


static int SelectDiscrete(int LTBF, int nDiscrGr, T_DiscrGr DiscrGr[mDiscrGr],
                          int mIx, int Ix[3])
{
  int    nIx, iIx;
  T_LTr  LLTr[mDiscrGr];
  int    nLLTr;

  if (nDiscrGr == 1) return 0;

  for (nIx = 1; nIx <= nDiscrGr - 1 && nIx <= mIx; nIx++)
  {
    range1(iIx, nIx) Ix[iIx] = iIx;

    do
    {
      ResetLLTr(LLTr, &nLLTr);
      range1(iIx, nIx)
        if (ExpLLTr(LTBF, mDiscrGr,
                    LLTr, &nLLTr, DiscrGr[Ix[iIx] + 1].P) < 0)
          return IE(-1);

      if (nLLTr >  nDiscrGr) return IE(-1);
      if (nLLTr == nDiscrGr) return nIx;
    }
    while (NextOf_n_from_m(nDiscrGr - 1, nIx, Ix) != 0);
  }

  return IE(-1);
}


static int CmpDiscr(const T_DiscrGr *a, const T_DiscrGr *b)
{
  return CmpiVect(a->Z, b->Z, 3);
}


static int Cmp_ssVM(const T_ssVM *a, const T_ssVM *b)
{
  return CmpiVect(a->V, b->V, 3);
}


int Set_ss(const T_SgOps *SgOps, T_ssVM ssVM[3])
{
  int        ir, i;
  int        nGen, IxGen[2];
  int        n_ssVM;
  T_RTMx     Z2PCBMx[2];
  int        nrSNF, DTBF, nd, id, d, f;
  int        SNF[3 * 3 * 3], Q[3 * 3], xp[3], x[3];
  int        nDiscrGr, iDG;
  T_LTr      LLTr[mDiscrGr];
  T_DiscrGr  DiscrGr[mDiscrGr];
  int        nIx, Ix[3];

  range1(ir, 3) rangei(3) ssVM[ir].V[i] = 0;
  range1(ir, 3) ssVM[ir].M = -1;

      nGen = SetAnyIxGen(SgOps, MGC_Undefined, IxGen);
  if (nGen < 0 || nGen > 2) return IE(-1);

      n_ssVM = GetContNullSpace(SgOps, IxGen, nGen, ssVM);
  if (n_ssVM < 0) return -1;
  if (n_ssVM == 3) return n_ssVM;

  if (GetZ2PCBMx(SgOps, Z2PCBMx) != 0) return -1;

      nrSNF = ConstructGenRmI(SgOps, Z2PCBMx, IxGen, nGen, SNF);
  if (nrSNF < 0) return IE(-1);

      nd = SmithNormalForm(SNF, nrSNF, 3, NULL, Q);
  if (nd < 0 || nd > 3) return IE(-1);

  DTBF = 1; range1(id, 3) DTBF = iLCM(DTBF, SNF[(nd + 1) * id]);

  ResetLLTr(LLTr, &nDiscrGr);

  range1(id, nd) {
    d = SNF[(nd + 1) * id];
    range2(f, 1, d) {
      rangei(3) xp[i] = 0;
      xp[id] = f * DTBF / d;
      iMxMultiply(x, xp, Q, 1, 3, 3);
      if (ExpLLTr(DTBF, mDiscrGr, LLTr, &nDiscrGr, x) < 0)
        return IE(-1);
    }
  }

  range1(iDG, nDiscrGr) {
    MemCpy(DiscrGr[iDG].P, LLTr[iDG].v, 3);
    RotMx_t_Vector(DiscrGr[iDG].Z, Z2PCBMx[1].s.R, DiscrGr[iDG].P, 0);
    rangei(3)
      DiscrGr[iDG].Z[i] = iModPositive(DiscrGr[iDG].Z[i], DTBF * CRBF);
  }

  if (BestVect(SgOps, ssVM, n_ssVM, DTBF, DiscrGr, nDiscrGr) != 0)
    return IE(-1);

  qsort((void *) DiscrGr, nDiscrGr, sizeof (*DiscrGr),
        (int (*)(const void *, const void *)) CmpDiscr);

      nIx = SelectDiscrete(DTBF, nDiscrGr, DiscrGr, 3 - n_ssVM, Ix);
  if (nIx < 0) return IE(-1);

  rangei(nIx) {
    if (n_ssVM >= 3) return IE(-1);
    MemCpy(ssVM[n_ssVM].V, DiscrGr[Ix[i] + 1].Z, 3);
    ssVM[n_ssVM].M = CancelBFGCD(ssVM[n_ssVM].V, 3, DTBF * CRBF);
    n_ssVM++;
  }

  qsort((void *) ssVM, n_ssVM, sizeof (*ssVM),
        (int (*)(const void *, const void *)) Cmp_ssVM);

  return n_ssVM;
}