#include <iostream.h>
#include "rules.h"
#include "rh_ilist.h"
#include "nn_bpg.h"
#include <math.h>

int    ms_get_iprop      (rx_eifc *eifc, const char* prop_name, int idx1=0, int idx2=0);
double ms_IP             (rx_eifc *eifc, int iels);
int    ms_lp_at_pis      (rx_eifc *eifc, int iel);
int    ms_deconjugate_pis(rx_eifc *eifc, int atom1, int atom2);


RULEOK ms_rules (ostream &print, rx_eifc *eifc, rx_datex *val, int irule, OP op, int sub_op, int test, int trace) {

static double     prob=0.;
static int const  *level=0, *huml=0;
static int        radic_els=0, dhuml=0;

static nn_bpg nn_alpha, nn_onium, nn_co;

switch(irule) {

case GLOBAL:
   switch(op) {
   case INIT_RULES: {
      val->putco("name","MS-rules for EROS7");
      val->putco("date","23.07.97");
      static const int       input_phase[]={ 0 };
      static const int phases_of_reactor[]={ 1, 0 };
      static const int    phase_property[]={ MONOMOLEC };
      static const int    phase_contacts[]={ 0 };
      static const int    input_together[]={ 0 };
      val->putco("input_phase_for_reactors", input_phase);
      val->putco("phases_per_reactor",       phases_of_reactor);
      val->putco("phase_property",           phase_property);
      val->putco("phase_contacts",           phase_contacts);
      val->putco("use_all_educts_together",  input_together);
      val->putco("output_phase",             1);
      static const char *const  kin_mode[]={ "prob_kin" };
      val->putco("kinetic_model",            kin_mode);
      val->putco("minimal_concentration",    1.e-5);
      val->put("reactivity",                 &prob);
      val->get("level",                      &level);
      val->get("huml",                       &huml);
      if (!huml) huml = &dhuml;
      if (!nn_alpha.init("rnnrul3x.dat")) return BAD;
      if (!nn_onium.init("rnnrul6zz.dat")) return BAD;
      if (!nn_co.init("rnnrul7z.dat")) return BAD;
      return OK;}
   case PREP_ROOT:
      return OK;
   case PREP_EDUCT: {
      radic_els = 0;
      int iel, nel, nelec;
      eifc->prop(nel, "E_N_ELECSYSS");
      for (iel = 1; iel<= nel; iel++) {
         eifc->prop(nelec, "EL_ELECTRONS",iel);
         if (nelec%2) {radic_els = iel; break;}
      }
      if (*level == 1) return OK;
      int charge=ms_get_iprop(eifc, "E_CHARGE");
      if (charge == 0) return BAD;
      return OK;}
   case DISTRIBFUNC:
      return OK;
   case FINISH:
      return OK;
   default:
      return BAD;
   }


case 1:
   switch(op) {
   case RULE_INFO: {
      static const char *const attrib[]={ "rearrangement", NULL };
      static const int         center[]={ 0 };
      val->putco("rule_name",            "Ionisation");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0: {
         if (ms_get_iprop(eifc, "E_CHARGE")) return BAD;
         return OK;}
      default:
         return BAD;
      }
   case FUNC: {
      int iel, nel, nelec, igr, sgr;
      rh_ilist<int> hgroup_idxs;
      vector<int>   hgroup;
      eifc->prop(nel, "E_N_ELECSYSS");
      for (iel = 1; iel <= nel; iel++) {
         eifc->set_active_ens(1);
         if (!ms_get_iprop(eifc, "EL_IS_SIGMA", iel)) {
            eifc->prop(hgroup, "EL_HGROUP", iel);
            igr = hgroup[0];
            sgr = hgroup[1];
            eifc->prop(nelec,"EL_ELECTRONS", iel);
            if ((!hgroup_idxs.at_idx(igr)) && (nelec)) {
               if (!ms_lp_at_pis(eifc, iel)) {
                  prob = 1e11 * exp(-2.*ms_IP(eifc, iel));
                  eifc->change_electrons(iel, -1);
                  eifc->symmetry_factor = sgr; // bzw. = sgr
print << "Ioni: prob = " << prob << endl;
                  eifc->save_active_ens_as_product();
                  hgroup_idxs.append(igr);
               }
            }
         }
      }
      return OK;}
   default:
      return BAD;
   }


case 2:
   switch(op) {
   case RULE_INFO: {
      static const char *const attrib[]={ NULL };
      static const int         center[]={ 1, 2, 0, 2, 3, 0, 0 };
      val->putco("rule_name",            "Alpha cleavage");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0: {
         if (!ms_get_iprop(eifc, "E_RADICAL")) return BAD;
         vector<int> atoms;
         if (!radic_els) return BAD;
         eifc->prop(atoms, "EL_ATOMS", radic_els);
         if (atoms.size() > 1) return BAD;
         return OK;}
      case 1:
         if (!ms_get_iprop(eifc, "A_RADICCENTER", eifc->center(1))) return BAD;
         return OK;
      case 2: {
         vector<int> els;
         eifc->prop(els, "N_ELECSYSS", eifc->center(1), eifc->center(2));
         if (els.size() > 2) return BAD;
         if ((els[0] == radic_els) || (els[els.size()] == radic_els)) return BAD;
         return OK;}
      case 3: {
         //if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(3)) == 1) return BAD;
         vector<int> els;
         eifc->prop(els, "N_ELECSYSS", eifc->center(2), eifc->center(3));
         if (els.size() > 2) return BAD;
         return OK;}
      default:
         return BAD;
      }
   case FUNC: {
      int issys, c1, c2, c3, iel1, iel2, imol, elem;
      double chrg;
      c1 = eifc->center(1);
      c2 = eifc->center(2);
      c3 = eifc->center(3);
      nn_alpha.set_vals(0, eifc);
      issys = ms_deconjugate_pis(eifc, c2, c3);
      if (!issys) return BAD;
      eifc->split_elsy(issys, c2, c3, iel1, iel2, 1);
      eifc->combine_elsys(radic_els, iel1, es_pi);
      eifc->prop(imol, "A_MOLECULE", c1);
      if (!ms_get_iprop(eifc, "M_CHARGE", imol)) return BAD;
      nn_alpha.set_vals(1, eifc);
      prob = nn_alpha.get_val();
      eifc->prop(elem, "A_ELEMENT", c3);
      eifc->prop(chrg, "A_CHARGECENTER", c1);
      if ((elem == 1) && (chrg < 0.1)) prob *= 0.06;
print << "alpha: prob = " << prob << endl;
      if (prob < 0.05) return BAD;
      return OK;}
   default:
      return BAD;
   }


case 3:
   switch(op) {
   case RULE_INFO: {
      static const char *const attrib[]={ NULL };
      static const int         center[]={ 1, 2, 0, 3, 4, 0, 0 };
      val->putco("rule_name",            "Onium reaction");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0:
         if (ms_get_iprop(eifc, "E_RADICAL")) return BAD;
         return OK;
      case 1: {
         int c1, elem, i, n, nel, nat, iel;
         double chrg;
         vector<int> els, atoms;
         c1 = eifc->center(1);
         eifc->prop(elem, "A_ELEMENT", c1);
         if ((elem == 1) || (elem == 6)) return BAD;
         eifc->prop(chrg, "A_CHARGECENTER", c1);
         if (chrg <= 0.) return BAD;
         eifc->prop(els, "A_ELECSYSS", c1);
         n = els.size();
         for (i=0; i<n; i++) {
            iel = els[i];
            if (!ms_get_iprop(eifc, "EL_IS_SIGMA", iel)) {
               eifc->prop(atoms, "EL_ATOMS", iel);
               nat = atoms.size();
               eifc->prop(nel, "EL_ELECTRONS", iel);
               if ((nat > 1) && (nel > 1)) return OK;
            }
         }
         return BAD;}
      case 2: {
         int c1=eifc->center(1), c2=eifc->center(2);
         if (ms_get_iprop(eifc, "A_ELEMENT", c2) != 6) return BAD;
         vector<int> els;
         eifc->prop(els, "N_ELECSYSS", c1, c2);
         if (els.size() != 1) return BAD;
         return OK;}
      case 3: {
         int c1=eifc->center(1), c2=eifc->center(2), c3=eifc->center(3), idist, idist2;
         if (ms_get_iprop(eifc, "A_ELEMENT", c3) != 6) return BAD;
         eifc->prop(idist, "N_DIST_INDEX", c2, c3);
         if (idist > 4) return BAD;
         eifc->prop(idist2, "N_DIST_INDEX", c1, c3);
         if (idist2 < idist) return BAD;
         return OK;}
      case 4:
         if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(4)) != 1) return BAD;
         return OK;
      default:
         return BAD;
      }
   case FUNC: {
      int c1, c2, c3, c4;
      c1 = eifc->center(1);
      c2 = eifc->center(2);
      c3 = eifc->center(3);
      c4 = eifc->center(4);
      nn_onium.set_vals(0, eifc);
      eifc->change_bond_order(c1, c2, -1);
      eifc->change_bond_order(c3, c4, -1);
      eifc->change_bond_order(c1, c4, 1);
      eifc->change_bond_order(c2, c3, 1);
      nn_onium.set_vals(1, eifc);
      prob = 0.5 * nn_onium.get_val();
print << "onium: prob = " << prob << endl;
      //if (prob < 0.05) return BAD;
      return OK;}
   default:
      return BAD;
   }


case 4:
   switch(op) {
   case RULE_INFO: {
      static const char *const attrib[]={ NULL };
      static const int         center[]={ 1, 2, 0, 2, 3, 0, 0 };
      val->putco("rule_name",            "Carbonyl elimination");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0:
         if (radic_els) return BAD;
         return OK;
      case 1: {
         double chrg;
         int elem = ms_get_iprop(eifc, "A_ELEMENT", eifc->center(1));
         if ((elem != 8) && (elem != 7)) return BAD;
         eifc->prop(chrg, "A_CHARGECENTER", eifc->center(1));
         if (chrg <= 0.) return BAD;
         return OK;}
      case 2: {
         int c1 = eifc->center(1), c2 = eifc->center(2);
         if (ms_get_iprop(eifc, "A_ELEMENT", c2) != 6) return BAD;
         vector<int> els;
         eifc->prop(els, "N_ELECSYSS", c1, c2);
         if (els.size() != 3) return BAD;
         return OK;}
      case 3:
         if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(3)) == 1) return BAD;
         return OK;
      default:
         return BAD;
      }
   case FUNC: {
      int issys, c2, c3, iel1, iel2;
      c2 = eifc->center(2);
      c3 = eifc->center(3);
      nn_co.set_vals(0, eifc);
      issys = ms_deconjugate_pis(eifc, c2, c3);
      if (!issys) return BAD;
      eifc->split_elsy(issys, c2, c3, iel1, iel2, 2);
      nn_co.set_vals(1, eifc);
      prob = 0.5 * nn_co.get_val();
print << "co: prob = " << prob << endl;
      //if (prob < 0.05) return BAD;
      return OK;}
   default:
      return BAD;
   }


case 5:
   switch(op) {
   case RULE_INFO: {
      static const char *const attrib[]={ NULL };
      static const int         center[]={ 1,2,0, 2,3,0, 3,4,0, 4,5,0, 5,6,0, 0 };
      val->putco("rule_name",            "McLafferty reaction");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0: {
         if (!radic_els) return BAD;
         vector<int> atoms;
         eifc->prop(atoms, "EL_ATOMS", radic_els);
         if (atoms.size() != 1) return BAD;
         return OK;}
      case 1: {
         int c1 = eifc->center(1), elem, i, n;
         eifc->prop(elem, "A_ELEMENT", c1);
         if (elem != 8) return BAD;
         vector<int> els;
         eifc->prop(els, "A_ELECSYSS", c1);
         n = els.size();
         for (i=0; i<n; i++) {
            if (els[i] == radic_els) return OK;
         }
         return BAD;}
      case 2: {
         vector<int> els;
         eifc->prop(els, "N_ELECSYSS", eifc->center(1), eifc->center(2));
         if (els.size() != 2) return BAD;
         return OK;}
      case 3:
         if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(3)) != 6) return BAD;
         return OK;
      case 4:
         return OK;
      case 5:
         return OK;
      case 6:
         if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(6)) != 1) return BAD;
         return OK;
      default:
         return BAD;
      }
   case FUNC: {
      int i, n, c1, c2, c3, c4, c5, c6, ipis, iel1, iel2, iel3, iel4, nhn, issys;
      vector<int> v;
      c1 = eifc->center(1);
      c2 = eifc->center(2);
      c3 = eifc->center(3);
      c4 = eifc->center(4);
      c5 = eifc->center(5);
      c6 = eifc->center(6);
      nhn = 0;
      eifc->prop(v, "A_NEIGHBORS", c5);
      n = v.size();
      for (i=0; i<n; i++) {
         if (ms_get_iprop(eifc, "A_ELEMENT", v[i]) != 1) nhn++;
      }
      issys = ms_deconjugate_pis(eifc, c3, c4);
      if (!issys) return BAD;
      eifc->prop(v, "N_ELECSYSS", c1, c2);
      n = v.size();
      for (i=0; i<n; i++) {
         if (!ms_get_iprop(eifc, "EL_IS_SIGMA", v[i])) ipis = v[i];
      }
      eifc->prop(v, "N_ELECSYSS", c5, c6);
      eifc->split_elsy(v[0], c5, c6, iel1, iel2, 1);
      eifc->split_elsy(issys, c3, c4, iel3, iel4, 1);
      eifc->combine_elsys(radic_els, iel2, es_sigma);
      eifc->combine_elsys(iel1, iel4, es_pi);
      eifc->combine_elsys(ipis, iel3, es_pi);
      prob = 0.8361;
      if (nhn == 1) prob = 0.23;
print << "ml: prob = " << prob << endl;
      //if (prob < 0.05) return BAD;
      return OK;}
   default:
      return BAD;
   }


case 6:
   switch(op) {
   case RULE_INFO: {
      if (!(*huml)) return BAD;
      static const char *const attrib[]={ "rearrangement", NULL };
      static const int         center[]={ 2, 3, 0, 0 };
      val->putco("rule_name",            "Hydrogen rearrangement");
      val->putco("attributes", attrib);
      val->putco("center_connectivity", center);
      return OK;}
   case CONSTR:
      switch(sub_op) {
      case 0: {
         if (!radic_els) return BAD;
         vector<int> atoms;
         eifc->prop(atoms, "EL_ATOMS", radic_els);
         if (atoms.size() != 1) return BAD;
         return OK;}
      case 1: {
         int c1 = eifc->center(1), i, n;
         vector<int> els;
         eifc->prop(els, "A_ELECSYSS", c1);
         n = els.size();
         for (i=0; i<n; i++) {
            if (els[i] == radic_els) return OK;
         }
         return BAD;}
      case 2: {
         int idist, c1 = eifc->center(1), c2 = eifc->center(2);
         eifc->prop(idist, "N_DIST_INDEX", c1, c2);
         if (ms_get_iprop(eifc, "A_ELEMENT", c1) == 6) {
            if ((idist < 3) || (idist > 5)) return BAD;
         } else {
            if (idist != 2) return BAD;
         }
         return OK;}
      case 3:
         if (ms_get_iprop(eifc, "A_ELEMENT", eifc->center(3)) != 1) return BAD;
         return OK;
      default:
         return BAD;
      }
   case FUNC: {
      int iel1, iel2, c1, c2, c3, elem, idist;
      vector<int> els;
      c1 = eifc->center(1);
      c2 = eifc->center(2);
      c3 = eifc->center(3);
      elem = ms_get_iprop(eifc, "A_ELEMENT", c1);
      eifc->prop(idist, "N_DIST_INDEX", c1, c2);
      eifc->prop(els, "N_ELECSYSS", c2, c3);
      eifc->split_elsy(els[0], c2, c3, iel1, iel2, 1);
      eifc->combine_elsys(radic_els, iel2, es_sigma);
      if (elem == 6) {
         switch(idist) {
         case 3:
            prob = 0.1;
            break;
         case 4:
            prob = 0.3;
            break;
         case 5:
            prob = 0.05;
            break;
         default:
            return BAD;
         }
      } else {
         prob = 0.5;
      }
      eifc->symmetry_factor = 1; // oder nicht?
print << "h-uml: prob = " << prob << endl;
      return OK;}
   default:
      return BAD;
   }


default:
   return BAD;


}

}


// --------------------------------------------------------------------------

int ms_get_iprop(rx_eifc *eifc, const char* prop_name, int idx1, int idx2) {
   int ivar=0;
   eifc->prop(ivar, prop_name, idx1, idx2);
   return ivar;
}

double ms_IP(rx_eifc *eifc, int iels) {
   return 10.;
}

int ms_lp_at_pis(rx_eifc *eifc, int iel) {
   int iel2, iat, nels, nat, i, j, nelec;
   if (ms_get_iprop(eifc, "EL_IS_SIGMA", iel)) return 0;
   vector<int> atoms, elsyss, atoms2;
   eifc->prop(atoms,"EL_ATOMS", iel);
   nat = atoms.size();
   if (nat == 1) return 0;
   for (i=0; i<nat; i++) {
      iat = atoms[i];
      eifc->prop(elsyss, "A_ELECSYSS", iat);
      nels = elsyss.size();
      for (j=0; j<nels; j++) {
         iel2 = elsyss[j];
         if (!ms_get_iprop(eifc,  "EL_IS_SIGMA", iel2)) {
            eifc->prop(atoms2, "EL_ATOMS", iel2);
            eifc->prop(nelec,  "EL_ELECTRONS", iel2);
            if ((nelec)&&(atoms2.size()==1)) return 1;
         }
      }
   }
   return 0;
}

int ms_deconjugate_pis(rx_eifc *eifc, int atom1, int atom2) {
   int issys=0, i, n, iel1, iel2;
   vector<int> elsys, etyp;
   eifc->prop(elsys, "N_ELECSYSS", atom1, atom2);
   n = elsys.size();
   etyp.resize(n);
   for (i=0; i<n; i++) {
      etyp[i] = ms_get_iprop(eifc, "EL_IS_SIGMA", elsys[i]);
   }
   for (i=0; i<n; i++) {
      if (etyp[i]) {
         issys = elsys[i];
      } else {
         if (eifc->split_elsy(elsys[i], atom1, atom2, iel1, iel2, -1)) return 0;
      }
   }
   return issys;
}