Cpi.py

#!/usr/bin/env python

# Copyright (C) 2019 Alireza Rashti #
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

from __future__ import division
import sympy
from sympy import *
from sympy.tensor.tensor import TensorIndexType, tensor_indices
from sympy.tensor.tensor import TensorHead, substitute_indices
from sympy import symbols, diag
from sympy import Eijk
from collections import OrderedDict
import re
import os

# import subprocess
import argparse
import datetime

# global vars:
CPI__glob_pr_flg = 0
CPI__glob_2dim_flg = 1
CPI__glob_Cpi_version = 3.0
CPI__glob_simplification_flag = -1
CPI__glob_index_stem = "___mu___"
__author__ = "Alireza Rashti"
__date__ = "June 2019"


def main():
  print("\nWelcome to Cpi\n")
  intro_print()

  # check the sympy version:
  version_digits = sympy.__version__.split(".")
  if float(version_digits[0]) < 1 or float(version_digits[1]) < 5:
    raise Exception(
        "\nThis script requires sympy version greater than or equal to 1.5.\n"
        "Your sympy version is {}".format(sympy.__version__))

  # read and parse inputs, e.g., file, flags, etc
  inputs = parse_inputs()

  # parsing system commands if any.
  sys_commands = parse_sys_commands(inputs)

  # generate C and python code
  file_name = gencode(inputs)

  # issue system commands:
  issue_system_commands(sys_commands, file_name)


# parsing system commands if any
def parse_sys_commands(arg):
  sys_commands = []

  for s in arg:
    if re.search(r"(?i)^Comm?and", s):
      s = re.sub(r'(?i)^Comm?and\["', "", s)
      s = re.sub(r'"\]@$', "", s)
      sys_commands.append(s)

  return sys_commands


# issue system commands on the output file:
def issue_system_commands(sys_commands, file_name):
  for c in sys_commands:
    cmd = c + " " + file_name
    print('Issuing command:\n"{}"\n'.format(cmd))
    os.system(cmd)

    # subprocess.call("ls")
    # subprocess.call([c,file_name])


# generating C code and python code according to the given input file
def gencode(arg):
  pr("#")
  print("Generating C code ...\n")
  # fill Maths info data base
  math_db = Maths_Info(arg)

  if CPI__glob_pr_flg:
    pr("-")
    math_db.pr()
    pr("-")

  # create C file with the name of input file
  C_file, file_name = creat_file(arg)

  # write and execute python code and realize C instructions due to calculations
  C_instructions = exec_pycode(math_db)

  # write the C instructions into the given file
  write_Ccode(C_instructions, math_db, C_file)

  # closing the file
  C_file.close()

  print("Generating C code --> done! :)\n")
  print('C code file is "{}".\n'.format(file_name))
  pr("#")

  return file_name


# write a C code according to the input file
def write_Ccode(C, CPI__db, Cfile):
  print("~~~~~~~> writing C code ...\n")

  # execute each C instructions:
  N = len(C)
  tab = ""
  now = datetime.datetime.now()
  ## NOTE: make sure the following copyright won't get longer than
  ##       4 lines including '/*' and newlines since some old cpi files
  ##       remove only these beginning 5 lines.
  intro = ("/*\n  These C codes generated by Cpi version {}\n"
           "  Copyright (C) 2019-{} Alireza Rashti.\n*/\n\n\n".format(
               CPI__glob_Cpi_version, now.year))
  fpr(Cfile, intro)

  # implement each instruction:
  for _ in range(N):
    # Ccode
    c = C[str(_)]
    if c["job"] == "Ccode":
      Cpr_str = c["Ccode"] + "\n"
      fpr(Cfile, Cpr_str)
      if re.search(r"^\\t+", c["Ccode"]):
        tab = re.search(r"^\\t+", c["Ccode"]).group(0)
      elif re.search(r"^\s+", c["Ccode"]):
        tab = re.search(r"^\s+", c["Ccode"]).group(0)
      else:
        tab = "  "

    # C declare
    elif c["job"] == "Cdeclare":
      declare_thingsC(CPI__db, Cfile, tab, c["declare_bin"])

    # populate
    elif c["job"] == "Cpopulate":
      Cpopulate(CPI__db, Cfile, c["Cpopulate"], tab)

    elif c["job"] == "Pcode":
      continue

    # calculations
    elif c["job"] == "calc":
      point = CPI__db.point_symb
      # if it has components
      if c["indexed"] == 1:
        sol = c[c["calc"]]
        for k in sol:
          ccode = "{0}double {1} =\n".format(tab, k)
          rhs = sol[k]

          # for the given fields that need to be evaluated on manifold point or has special arguments (C_arg)
          for symb in CPI__db.symbols_ld:
            if symb["obj"] == "field":
              if "C_argument" in symb.keys():
                array = set(symb["array_comp"]) # remove the redundants
                for e in array:
                  if e != "0." and not re.search(r"^-", e):
                    Carg = symb[symb["C_argument"]]
                    Carg = re.sub(r"name", e, Carg)
                    ################ if there is U? or D? ###################
                    if re.search(r"(?i)[UD]+\?", Carg):
                      Carg_indices = re.findall(r"[UD]+\?", Carg)
                      lCarg_indices = len(Carg_indices)
                      # substitute the question mark:
                      for _ in range(lCarg_indices):
                        Carg_indices[_] = re.sub(r"\?", "", Carg_indices[_])
                      # some checks rank and type:
                      if lCarg_indices != int(symb["rank"]): # check rank agrees
                        raise Exception(
                            "The indices for C_arg are not matched between {} and {}.\n"
                            .format(
                                symb["name"],
                                symb[symb["C_argument"]],
                            ))
                      t0 = 0
                      for t in symb["type"]: # check type agrees
                        if not re.search(
                            r"(?i){}".format(Carg_indices[t0]),
                            t,
                        ):
                          raise Exception(
                              "The indices for C_arg are not matched between {} and {}.\n"
                              .format(
                                  symb["name"],
                                  symb[symb["C_argument"]],
                              ))
                        else:
                          t0 += 1
                      # parse the indices of component:
                      suffix = re.search(r"_[UD\d]+$", e).group(0)
                      # suffix = re.sub(r'^_','',suffix)
                      indices = re.findall(r"\d", suffix)
                      # now substitue indices in Carg
                      Carg_split = Carg.split(",")
                      lCarg_split = len(Carg_split)
                      c1 = 0
                      c2 = 0
                      Carg2 = ""
                      for p in Carg_split:
                        if re.search(r"(?i)[UD]+\?", p):
                          Carg_split[c1] = re.sub(
                              r"(?i)[UD]+\?",
                              "{}".format(indices[c2]),
                              p,
                          )
                          c2 += 1
                        if c1 < lCarg_split - 1:
                          Carg2 += (Carg_split[c1] + ",")
                        c1 += 1

                      Carg2 += Carg_split[c1 - 1]
                      Carg = Carg2
                    ################ END of if there is U? or D? ###################
                    rhs = re.sub(
                        r"\b{}\b".format(e),
                        "{}{}".format(e, Carg),
                        rhs,
                    )

              else: # if there is no C_arg
                array = set(symb["array_comp"]) # remove the redundants
                for e in array:
                  if e != "0." and not re.search(r"^-", e):
                    rhs = re.sub(
                        r"\b{}\b".format(e),
                        "{}[{}]".format(e, point),
                        rhs,
                    )

          ccode += "{};\n".format(rhs)
          ccode = re.sub(r"\b1.0\*\b", "", ccode)
          ccode = style_eq_pr(ccode, tab)

          fpr(Cfile, ccode)

      # if doesn't have components
      else:
        ccode = "{0}double {1} =\n".format(tab, c["calc"])
        rhs = c[c["calc"]]

        # for the given fields that need to be evaluated on manifold point or has special arguments (C_arg)
        for symb in CPI__db.symbols_ld:
          if symb["obj"] == "field":
            if "C_argument" in symb.keys():
              array = set(symb["array_comp"]) # remove the redundants
              for e in array:
                if e != "0." and not re.search(r"^-", e):
                  Carg = symb[symb["C_argument"]]
                  Carg = re.sub(r"name", e, Carg)
                  ################ if there is U? or D? ###################
                  if re.search(r"(?i)[UD]+\?", Carg):
                    Carg_indices = re.findall(r"[UD]+\?", Carg)
                    lCarg_indices = len(Carg_indices)
                    # substitute the question mark:
                    for _ in range(lCarg_indices):
                      Carg_indices[_] = re.sub(r"\?", "", Carg_indices[_])
                    # some checks rank and type:
                    if lCarg_indices != int(symb["rank"]): # check rank agrees
                      raise Exception(
                          "The indices for C_arg are not matched between {} and {}.\n"
                          .format(
                              symb["name"],
                              symb[symb["C_argument"]],
                          ))
                    t0 = 0
                    for t in symb["type"]: # check type agrees
                      if not re.search(
                          r"(?i){}".format(Carg_indices[t0]),
                          t,
                      ):
                        raise Exception(
                            "The indices for C_arg are not matched between {} and {}.\n"
                            .format(
                                symb["name"],
                                symb[symb["C_argument"]],
                            ))
                      else:
                        t0 += 1
                    # parse the indices of component:
                    suffix = re.search(r"_[UD\d]+$", e).group(0)
                    # suffix = re.sub(r'^_','',suffix)
                    indices = re.findall(r"\d", suffix)
                    # now substitue indices in Carg
                    Carg_split = Carg.split(",")
                    lCarg_split = len(Carg_split)
                    c1 = 0
                    c2 = 0
                    Carg2 = ""
                    for p in Carg_split:
                      if re.search(r"(?i)[UD]+\?", p):
                        Carg_split[c1] = re.sub(
                            r"(?i)[UD]+\?",
                            "{}".format(indices[c2]),
                            p,
                        )
                        c2 += 1
                      if c1 < lCarg_split - 1:
                        Carg2 += Carg_split[c1] + ","
                      c1 += 1

                    Carg2 += Carg_split[c1 - 1]
                    Carg = Carg2
                  ################ END of if there is U? or D? ###################
                  rhs = re.sub(
                      r"\b{}\b".format(e),
                      "{}{}".format(e, Carg),
                      rhs,
                  )

            else: # if there is no C_arg
              array = set(symb["array_comp"]) # remove the redundants
              for e in array:
                if e != "0." and not re.search(r"^-", e):
                  rhs = re.sub(
                      r"\b{}\b".format(e),
                      "{}[{}]".format(e, point),
                      rhs,
                  )

        ccode += "{};\n".format(rhs)
        ccode = re.sub(r"\b1.0\*\b", "", ccode)
        ccode = style_eq_pr(ccode, tab)

        fpr(Cfile, ccode)

    else:
      raise Exception("No job!")


# populating component in C
def Cpopulate(CPI__db, C_file, pop, tab):
  fpr(C_file, "\n")
  sp = pop.split("=")
  lhs = sp[0]
  rhs = sp[1]

  # check both lhs have same rank and symmetry
  lhs_obj = []
  rhs_obj = []
  flg1 = 1
  flg2 = 1

  for obj in CPI__db.symbols_ld:
    if obj["name"] == lhs and flg1:
      lhs_obj = obj
      flg1 = 0
    elif obj["name"] == rhs and flg2:
      rhs_obj = obj
      flg2 = 0
    elif not flg1 and not flg2:
      break

  try:
    if rhs_obj["obj"] != "local":
      raise Exception(
          "The Right hand side in Cpopulate[{}] cannot be used for assignment.".
          format(pop))
  except:
    raise Exception(
        "The Right hand side in Cpopulate[{}] cannot be used for assignment.".
        format(pop))

  if lhs_obj["rank"] != rhs_obj["rank"]:
    raise Exception("Ranks are not equal for Cpopulate[{}].".format(pop))

  if "symmetry" in rhs_obj.keys() or "symmetry" in lhs_obj.keys():
    try:
      # check all symmetries exist
      for s in lhs_obj["symmetry"]:
        if s not in rhs_obj["symmetry"]:
          raise Exception(
              "Symmetries are not equal for Cpopulate[{}].".format(pop))
      for s in rhs_obj["symmetry"]:
        if s not in lhs_obj["symmetry"]:
          raise Exception(
              "Symmetries are not equal for Cpopulate[{}].".format(pop))
    except:
      raise Exception("Symmetries are not equal for Cpopulate[{}].".format(pop))

  # now populate:
  lhs_cmp = ordered(set(lhs_obj["array_comp"]))

  for cmp in lhs_cmp:
    if cmp != "0." and not re.search(r"^-", cmp):
      suffix = re.sub(r"^{}".format(lhs_obj["name"]), "", cmp)
      rhs_cmp = rhs_obj["name"] + suffix
      ccode = (tab + cmp + "[" + CPI__db.point_symb + "]" + " = " + rhs_cmp +
               ";" + "\n")
      fpr(C_file, ccode)


# write and execute the python code according to the input file
# and realize C instructions due to calculations
# NOTE: exec() function has access to the local variable of the
# exec_pycode() function. as a result, it might happen that this
# variables have the similar name as what written in the input file!
# to fix this, each variable's name should have peculiar name.
def exec_pycode(CPI__db):
  print("~~~~~~~> writing and executing sympy code ...\n")

  py_head = ""
  C_instructions = dict()
  repl = "Tensors_db"
  # declare symbols
  for s in CPI__db.symbols_ld:
    if ("KD" == s["name"] or "EIJK" == s["name"] or
        re.search(r"KD.+i", s["name"])):
      continue

    if "array_comp" in s.keys():
      arry_set = sorted(set(s["array_comp"])) # reduce the redundants

      for comp in arry_set:
        # No 0 and - sign comming from symmetry
        if comp != "0." and not re.search(r"^-", comp):
          py_head += "{0:20} = symbols('{0}')\n".format(comp)

  # declare functions:
  for s in CPI__db.symbols_ld:
    if s["obj"] == "function":
      py_head += "{0:20} = Function('{0}')\n".format(s["name"])

  # declare indices
  py_head += 'L = TensorIndexType("L",dim = {})\n'.format(CPI__db.dim_i)
  for ind in CPI__db.indices_s:
    py_head += '{0:5} = tensor_indices("{0}",L)\n'.format(ind)

  # declare tensors
  for obj in CPI__db.symbols_ld:
    if (obj["obj"] == "field" or
        obj["obj"] == "local") and obj["type"] != "scalar":
      n = int(obj["rank"])
      L = "["
      ones = "["
      for i in range(n - 1):
        L += "L,"
        ones += "[1],"
      L += "L]"
      ones += "[1]]"
      #      py_head +='{0:20} = TensorHead("{0}",{1},{2})\n'.format(obj['name'],L,ones)
      #     the above line was used for sympy 1.4
      py_head += '{0:20} = TensorHead("{0}",{1})\n'.format(obj["name"], L)

  # populate tensors
  L = "diag("
  for i in range(CPI__db.dim_i - 1):
    L += "1.,"
  L += "1.)"
  indices = indices_tuple(CPI__db, CPI__db.dim_i)

  py_head += "{} = {{L:{}}}\n".format(repl, L)

  for obj in CPI__db.symbols_ld:
    if (obj["obj"] == "field" or
        obj["obj"] == "local") and obj["type"] != "scalar":
      matrix = "{}".format(obj["matrix_comp"])
      matrix = re.sub(r"'", "", matrix)
      indices = indices_tuple(CPI__db, int(obj["rank"]))
      py_head += "{0}.update({{{1}{2}:{3}}})\n".format(repl,
                                                       obj["name"],
                                                       indices,
                                                       matrix)

  n = len(CPI__db.instruction_dd)
  C_instructions = dict() # C instructions
  for _i in range(n):
    py_code = py_head

    # evaluate equations and expand indices
    eq_sol = ""
    extra_ind = set()

    if CPI__db.instruction_dd[str(_i)]["job"] == "Ccode":
      job = dict()
      job["job"] = "Ccode"
      job["Ccode"] = CPI__db.instruction_dd[str(_i)]["Ccode"]
      C_instructions[str(_i)] = job

    elif CPI__db.instruction_dd[str(_i)]["job"] == "Cdeclare":
      job = dict()
      job["job"] = "Cdeclare"
      job["declare_bin"] = CPI__db.instruction_dd[str(_i)]["declare_bin"]
      C_instructions[str(_i)] = job

    elif CPI__db.instruction_dd[str(_i)]["job"] == "Cpopulate":
      job = dict()
      job["job"] = "Cpopulate"
      job["Cpopulate"] = CPI__db.instruction_dd[str(_i)]["Cpopulate"]
      C_instructions[str(_i)] = job

    elif CPI__db.instruction_dd[str(_i)]["job"] == "Pcode":
      job = dict()
      job["job"] = "Pcode"
      C_instructions[str(_i)] = job
      py_head += CPI__db.instruction_dd[str(_i)]["Pcode"]
      py_head += "\n"

    elif CPI__db.instruction_dd[str(_i)]["job"] == "calc":
      job = dict()
      job["job"] = "calc"
      calc = CPI__db.instruction_dd[str(_i)]
      eq = calc["calc"].split("=")
      lhsO = eq[0] # O for original
      rhsO = eq[1] # O for original
      lhsA = lhsO # A for Amenable
      rhsA = rhsO # A for Amenable

      print("-> eq: {}  ...\n".format(calc["calc"]))

      # if lhs is tensorial:
      if re.search(r"\([-\w,]+\)", lhsA):
        lhsA = re.sub(r"\([-\w,]+\)", "", lhsA)

      eqA = "{:14}".format(lhsA) + "= " + rhsA + "\n"

      # making C instructions the syntax is as follows:
      # C_instructions[job_number][lhs_name][lhs_component] = rhs_components
      job["calc"] = lhsA
      subjob = dict()
      arry_set = set()
      lhs_obj = dict()
      comp_flg = 0
      lhs_rank = 0
      for s in CPI__db.symbols_ld:
        if s["name"] == lhsA and "array_comp" in s.keys():
          arry_set = set(s["array_comp"]) # reduce the redundants
          lhs_obj = s
          lhs_rank = int(s["rank"])
          comp_flg = 1

      if comp_flg == 1:
        for comp in arry_set:
          # 0 and - sign comming from symmetry
          if not re.search(r"^0", comp) and not re.search(r"^-", comp):
            subjob[comp] = ""

      # NOTE: in python 2 the order is NOT kept when new entries
      # are added to a dict.
      subjob = OrderedDict(sorted(subjob.items()))

      sol = "" # C solution
      unmatched_indices = list()
      # if lhs is not tensorial
      if not re.search(r"\([-\w,]+\)", lhsO):
        job["indexed"] = 0
        sol = 'try:\n\tC_instructions["{0}"]["{1}"] = ccode(my_simplify(({1}.substitute_indices().replace_with_arrays({2})).doit()))\n'.format(
            str(_i), lhsA, repl)
        #        sol += 'except 1:\n\tC_instructions["{0}"]["{1}"] = ccode(my_simplify({1}.replace_with_arrays({2})))\n'.format(str(_i),lhsA,repl)
        #        sol += 'except AttributeError:\n\tC_instructions["{0}"]["{1}"] = ccode(my_simplify({1}.doit()))\n'.format(str(_i),lhsA,repl)
        sol += 'except:\n\tC_instructions["{0}"]["{1}"] = ccode(my_simplify({1}))\n'.format(
            str(_i), lhsA, repl)
        #        sol += 'except:\n\tC_instructions["{0}"]["{1}"] = ccode(my_simplify({1}))\n'.format(str(_i),lhsA,repl)
        sol += "try:\n\tres = str(type({}.replace_with_arrays({})))\n".format(
            lhsA, repl)
        sol += "\tif 'ImmutableDenseNDimArray' in res:\n"
        sol += "\t\tunmatched_indices.append('yes')\n"
        sol += "except:\n\tunmatched_indices.append('no')\n"

      # if lhs is tensorial
      else:
        assert comp_flg == 1

        job["indexed"] = 1
        ## get rhs free indices, ex: g(i,j) = ... => repl_indices = [i,j]
        repl_indices = re.search(r"\([-\w,]+\)",
                                 lhsO).group(0) # only - sign allowed
        repl_indices = re.sub(r"^\(", "[", repl_indices)
        repl_indices = re.sub(r"\)$", "]", repl_indices)
        repl_indices = re.sub(r"-", "", repl_indices)

        for comp in arry_set:
          # note 0 and - sign comming from symmetry
          if not re.search(r"^0", comp) and not re.search(r"^-", comp):
            # getting indices:
            e = re.sub(r"{}_".format(lhsA), "", comp)
            e = re.sub(r"[UD]", " ", e)
            index = e.split(" ")
            index.remove("")
            n = len(index)
            suffix = "{0}.substitute_indices().replace_with_arrays({1},{2})[".format(
                lhsA, repl, repl_indices)
            for i in range(n - 1):
              suffix += "{},".format(int(index[i]))
            suffix += "{}]".format(int(index[n - 1]))

            sol += 'try:\n\tC_instructions["{0}"]["{1}"]["{2}"] = ccode(my_simplify(({3}).doit()))\n'.format(
                str(_i), lhsA, comp, suffix)
            sol += 'except:\n\tC_instructions["{0}"]["{1}"]["{2}"] = ccode(my_simplify({3}))\n'.format(
                str(_i), lhsA, comp, suffix)

      eqA = "{:14}".format(lhsA) + "= " + rhsA + "\n"
      eq_sol += eqA + sol

      job[lhsA] = subjob
      C_instructions[str(_i)] = job

      for ind in extra_ind:
        py_code += '{0:5} = tensor_indices("{0}",L)\n'.format(ind)

      py_code += eq_sol

      if CPI__glob_pr_flg:
        pr("-")
        print("sympy code for {}:".format(calc["calc"]))
        print(py_code)
        pr("-")

      exec(py_code)

      if len(unmatched_indices) == 1:
        if unmatched_indices[0] == "yes":
          raise Exception(
              "\nIn eq: {}\nthe free indices are not matched!\n".format(
                  calc["calc"]))

      try:
        print("-> eq: {}    {}\n".format(calc["calc"], "\u2713"))
      except:
        print("-> eq: {}  done.\n".format(calc["calc"]))

    else:
      raise Exception("No job")

  return C_instructions


# declare thins in C file as it is given from the input
def declare_thingsC(CPI__db, C_file, tab, bin):
  # fpr(C_file,'\n')

  for obj in CPI__db.symbols_ld:
    # for each user call declare we have declare_bin thus:
    if not "declare_bin" in obj.keys():
      continue
    # if bin is not matched skip
    if obj["declare_bin"] != bin:
      continue

    # variables:
    if obj["obj"] == "variable":
      if obj["Ccall"] == "Ccode":
        fpr(C_file, obj["Ccode"])
      elif re.search(r"C_macro", obj["Ccall"]):
        Cstr = C_macro_replace_built_in(obj[obj["Ccall"]],
                                        obj["name"],
                                        obj["name"],
                                        "0",
                                        tab)
        fpr(C_file, Cstr)

      elif obj["Ccall"] == "none":
        continue
      else:
        raise Exception("\n No method to declare {} in C.".format(obj["name"]))

    # scalars:
    elif obj["obj"] == "field" and obj["type"] == "scalar":
      if obj["Ccall"] == "Ccode":
        fpr(C_file, obj["Ccode"])
      elif re.search(r"C_macro", obj["Ccall"]):
        Cstr = C_macro_replace_built_in(obj[obj["Ccall"]],
                                        obj["name"],
                                        obj["name"],
                                        "0",
                                        tab)
        fpr(C_file, Cstr)

      elif obj["Ccall"] == "none":
        continue
      else:
        raise Exception("\n No method to declare {} in C.".format(obj["name"]))

    # tensors:
    elif obj["obj"] == "field" and obj["type"] != "scalar":
      if not "Ccall" in obj.keys():
        raise Exception("\n No method to declare {} in C.".format(obj["name"]))
      elif obj["Ccall"] == "Ccode":
        fpr(C_file, obj["Ccode"])
      elif re.search(r"C_macro", obj["Ccall"]):
        array = sorted(set(obj["array_comp"]))
        comp_counter = 0
        for cmp in array: # for each component
          if cmp != "0." and not re.search(r"^-", cmp):
            Cstr = C_macro_replace_built_in(
                obj[obj["Ccall"]],
                cmp,
                obj["name"],
                str(comp_counter),
                tab,
            )
            fpr(C_file, Cstr)

            comp_counter += 1

      elif obj["Ccall"] == "none":
        continue
      else:
        raise Exception("\n No method to declare {} in C.".format(obj["name"]))

  # fpr(C_file,"\n\n")


# populating components of the given objects in matrix and array
def populate_components(obj, CPI__db):
  if obj["obj"] == "variable":
    return [obj["name"]], [obj["name"]]

  elif (obj["obj"] == "field" or
        obj["obj"] == "local") and obj["type"] == "scalar":
    return [obj["name"]], [obj["name"]]

  elif (obj["obj"] == "field" or
        obj["obj"] == "local") and obj["type"] != "scalar":
    return realize_components(obj, CPI__db)
  else:
    return [], []


# given an object it returns list of components of the object according
# to the symmetry of the objects and type and rank in array format and indexed obj format
def realize_components(obj, CPI__db):
  name = obj["name"]

  if not "rank" in obj.keys():
    raise Exception("This object {} does not have rank.\n".format(obj["name"]))

  if not "type" in obj.keys():
    raise Exception("This object {} does not have type.\n".format(obj["name"]))

  if re.search(r"^\d", name):
    raise Exception("{} should not start with a number.\n".format(name))

  rank = int(obj["rank"])
  type = obj["type"]
  tab = ""
  comp = "'" + name + "_"
  append = "'.format("
  append2 = "["
  indx = 0
  code = ""

  for i in range(rank): # fors
    if i > 0:
      append += ","
      append2 += ","
      tab += "\t"
    comp += "{}{{}}".format(type[i])
    append += "arg[{}]".format(i)
    append2 += "_{}".format(i)
    if i > 0:
      code += "{}comp{} = []\n".format(tab, i)
    code += "{}for _{} in range({}):\n".format(tab, i, CPI__db.dim_i)
    indx = i

  # populate the whole tensor with no symmetry
  tab += "\t"
  append += ")"
  append2 += "]"
  code += "{}arg = {}\n".format(tab, append2)
  code += "{}comp{} = {}{}\n".format(tab, indx + 1, comp, append)
  code += "{}array.append(comp{})\n".format(tab, indx + 1)
  code += "{}comp{}.append(comp{})\n".format(tab, indx, indx + 1)
  # going upward of the nested loop
  for i in range(rank - 1, 0, -1):
    tab = re.sub(r"\s{1}?$", "", tab)
    code += "{}comp{}.append(comp{})\n".format(tab, i - 1, i)
  comp0 = [] # matrix format of the indexed_obj (misnomer)
  array = [] # array format of the the indexed_obj
  try:
    exec(code)
  except:
    raise Exception("Cannot populate components of '{}'!".format(name))

  # if there is some symmetry
  if len(obj["symmetry"]) != 0:
    # loop over all symmetries and modify the components accordingly
    for s in obj["symmetry"]:
      comp0, array = reduce_symm_components(comp0, s, obj, CPI__db)

  indexed_obj = comp0
  return indexed_obj, array


# realize the given equation in python
def realize_eq_py(CPI__db, eq):
  sympy_defined = [] # dir(sympy)
  sympy_defined.append("L")
  # checkup for symbols and add new symbols if necessary
  parts = eq.split("=")
  lhs = parts[0]
  rhs = parts[1]

  # check if the symbols in lhs is not already defined
  name = re.sub(r"(?<=\w)\({1}(-?[a-zA-Z]+\w?,?)+\){1}", "", lhs)
  for obj in CPI__db.symbols_ld:
    if obj["name"] == name:
      print("{} has been defined more than one time.\n"
            "For each new definition one must use new name.\n".format(lhs))
      raise Exception("{}={}\n{} is protected.\n".format(lhs, rhs, obj["name"]))

  if re.search(r"[\.\+\*/]", lhs):
    raise Exception(
        "It cannot realize the left hand side of eq:\n{}.".format(eq))

  # check if there is rank 3 or higher tensor on rhs with lower indices(covariant)
  # as I observed there is a bug in sympy 1.5dev version which cannot correctly
  # convert upper indices (contravariant) tensor of rank higher than 2 to lower indices (covariant).
  # if they solve this issue, this check can be removed. right now, the way we deal with higher
  # rank tensor of 3 or more would be with Kronecker delta (KD). So in the calculation
  # we always use contravariant indices and we contract them with using KD as the metric.
  for obj in CPI__db.symbols_ld:
    if "rank" in obj.keys():
      if int(obj["rank"]) > 2:
        if re.search(
            r"\b{}\({{1}}[\w,\-]*\-[\w,\-]*\){{1}}".format(obj["name"]),
            rhs,
        ):
          raise Exception(
              "\nEq: {} got problem.\nPlease read the following:\n".format(eq) +
              "I observed there is a bug in sympy 1.4 and higher versions which cannot correctly\n"
              "convert upper indices (contravariant) tensor of rank higher than 2 to lower indices (covariant).\n"
              "If they solve this issue, this check can be removed. right now, the way we deal with higher\n"
              "rank tensor of 3 or more will be with Kronecker delta (KD). So in the calculation\n"
              "we always use contravariant indices and we contract them with using KD as the metric.\n"
              "For example:\n"
              "cont = x(i,j,k)*x(-i,-j,-k) -> cont = x(i,j,k)*x(l,m,n)*KD(-i,-l)*KD(-j,-m)*KD(-k,-n).\n"
          )

  # now since this object is not defined let's add it to the CPI__db
  new_obj = dict()
  new_obj["name"] = name
  if name in sympy_defined:
    raise Exception(
        "The name {} is protected for Sympy library, use another name.".format(
            name))

  new_obj["obj"] = "local"
  rank = 0
  if re.search(r"\(.*\)", lhs):
    rank = len(lhs.split(",")) # e.g psi(i,j) would be 2
  new_obj["rank"] = str(rank)
  # type:
  if rank != 0:
    ty = re.sub(r"^\w+\(", "", lhs)
    ty = re.sub(r"\)$", "", ty)
    piece = ty.split(",")
    type = []
    for p in piece:
      if re.search(r"^-\w+", p):
        type.append("D")
      else:
        type.append("U")

    new_obj["type"] = type
  else:
    new_obj["type"] = "scalar"

  CPI__db.symbols_ld.append(new_obj)

  # get all of the symbols in rhs by cutting their indices
  trimmed = rhs
  # exception is for integrate, which might need extra symbol
  # so I cut it off completly
  if "integrate" in rhs:
    trimmed = re.sub(r"integrate\([\w\-\*\+/\(\)]+,\(?[\w\-\*\+/,]+\)?\)",
                     "1.",
                     trimmed)

  # get rid of tensor indices
  trimmed = re.sub(r"(?<=\w)\({1}(-?[a-zA-Z]+\w?,?)+\){1}", "", trimmed)
  # get rid of sympy functions and pars their argumnets
  for s in sympy_defined:
    if re.search(r"\b{}\b".format(s), rhs):
      # if function:
      if re.search(r"\b{}\(".format(s), trimmed):
        trimmed = re.sub(r"\b{}".format(s), "1.*", trimmed)
        trimmed = re.sub(r",".format(s), "+", trimmed)
      # if constant number:
      if re.search(r"\b{}\b".format(s), trimmed):
        trimmed = re.sub(r"{}".format(s), "1.", trimmed)

  # some light check if the symbols in rhs all defined
  # there might some quantites scape from this check but will be catch later
  code = str()
  N = len(CPI__db.symbols_ld)
  for _ in range(N - 1):
    obj = CPI__db.symbols_ld[_]
    code += '{0:10} = symbols("{0}")\n'.format(obj["name"])

  # add kronecker delta and levi civita too
  code += '{0:10} = symbols("{0}")\n'.format("KD")
  code += '{0:10} = symbols("{0}")\n'.format("KD0i")
  code += '{0:10} = symbols("{0}")\n'.format("KD1i")
  code += '{0:10} = symbols("{0}")\n'.format("KD2i")
  code += '{0:10} = symbols("{0}")\n'.format("KD3i")
  code += '{0:10} = symbols("{0}")\n'.format("EIJK")
  code += "expr = {}\n".format(trimmed)
  code += "___t___  = srepr(expr)\n"

  if trimmed != "":
    try:
      exec(code)
    except:
      raise Exception(
          "An error occurred in\n{}."
          " It might be for either undefined variables or syntax error.".format(
              eq))


# given file and string it prints the str into the file
def fpr(f, str):
  # take care of \n
  if re.search(r".*\\n+.*", str):
    if re.search(r'(?:").*\\n.*(?:")', str):
      str_quot = re.search(r'(?:").*\\n.*(?:")', str).group(0)
      split_quot = str.split(str_quot)
      if len(split_quot) == 2:
        s0 = split_quot[0].split("\\n")
        s0_new = "\n".join(s0)
        s1 = split_quot[1].split("\\n")
        s1_new = "\n".join(s1)
        str = s0_new + str_quot + s1_new
      else:
        raise Exception(
            "It cannot realize \\n in {} due to quotation mark.\n".format(str))
    else:
      str2 = str.split("\\n")
      str = "\n".join(str2)
  # take care of \t
  if re.search(r".*\\t+.*", str):
    str2 = str.split("\\t")
    str = "  ".join(str2)

  f.write(str)


# given a symmetry relation, it finds the kind of symmetry
# and gives a dictionary after realization
def realize_symmetry(relation):
  # check the format:
  if not re.search(r"^[-\w\(,\)]+=[-+]?[-\w\(,\)]+$", relation):
    raise Exception(
        "The form of symmetry relation {} is not correct.".format(relation))

  # symm or anti-symm
  sign = ""
  if re.search(r"=-{1}", relation):
    sign = "-"
  else:
    sign = "+"

  # find places of permuted indices
  perm = []
  sym = dict()
  l1 = dict()
  l2 = dict()
  lhs = re.search(r"^[-\w\(,\)]+[^=]", relation).group(0)
  rhs = re.search(r"(?:=)[-+]?[-\w\(,\)]+", relation).group(0)
  rhs = re.sub(r"=", "", rhs)

  # test
  # print(lhs)
  # print(rhs)
  # end

  # trims
  lhs = re.sub(r"^\w+", "", lhs)
  rhs = re.sub(r"^[-+]?\w+", "", rhs)
  lhs = re.sub(r"[\(\)]", "", lhs)
  rhs = re.sub(r"[\(\)]", "", rhs)
  lhs_l = lhs.split(",")
  rhs_l = rhs.split(",")

  # test
  # print(lhs_l)
  # print(rhs_l)
  # end

  # checkup
  if len(lhs_l) != len(rhs_l):
    raise Exception("The indices in {} are not matches!".format(relation))

  for _0 in lhs_l:
    checkup_flg = 0
    for _1 in rhs_l:
      if _0 == _1:
        checkup_flg = 1
        break
    if checkup_flg == 0:
      raise Exception("The indices in {} are not matches!".format(relation))

  n = len(lhs_l)
  for i in range(n):
    if rhs_l[i] != lhs_l[i]:
      perm.append(i)

  # if the symmetry relation is trivial like g(i,j):= g(i,j)
  if len(perm) == 0:
    raise Exception("No two indices are changed in {}.".format(relation))
  # checkup
  if len(perm) > 2:
    raise Exception(
        "Only two indices can be permuted for symmetry in {}!".format(relation))

  # add sign info
  sym["sign"] = sign
  sym["perm"] = perm

  # test
  # print(sym)
  # end

  return sym


# print custom line:
def pr(s):
  str = ""
  for i in range(70):
    str += s
  str += "\n"

  print(str)


# creating files
def creat_file(arg):
  flg = 0
  for s in arg:
    if re.search(r"^file_name=\w+\..+$", s):
      file_name = re.sub(r"^file_name=", "", s)
      file_name = re.sub(r"\..+$", "", file_name)
      flg += 1
  if flg == 0:
    raise Exception("It could not find the file name!")

  # making a foldor for calculations
  path = os.getcwd()
  i = 0
  dir = path + "/{}.c".format(file_name, i)
  while os.path.exists(dir):
    i += 1
    dir = path + "/{}_{:02}.c".format(file_name, i)

  if i > 0:
    file_name = "{}_{:02}".format(file_name, i)

  assert file_name

  file_name += ".c"
  C_file = open(file_name, "wt")

  return C_file, file_name


# check up the input if there is any error
def checkup_input(arg):
  return
  # check for redundancy
  n = len(arg)
  for s in arg:
    for i in range(arg.index(s) + 1, n):
      if s == arg[i] and ("Ccode" not in arg[i] and "ccode" not in arg[i]):
        raise Exception(
            'Duplicated parameter "{}" in the input file!'.format(s))


# remove white spaces and comments
def trim_inputs(arg):
  arg_new = []

  # join broken lines
  N = len(arg)
  i = 0
  while i < N:
    if re.search(r"\\\n$", arg[i]):
      m = ""
      while re.search(r"\\\n$", arg[i]) and i < N - 1:
        arg[i] = re.sub(r"\\\n$", "", arg[i])
        m += arg[i]
        i += 1
      m += arg[i]
      arg_new.append(m)
    else:
      arg_new.append(arg[i])
    i += 1

  # remove newlines and comments
  arg = arg_new
  pat_newline = re.compile(r"^\s*\n$")
  pat_comment = re.compile(r"^\s*#")
  n = len(arg)
  p = 0
  while p < n:
    if pat_newline.search(arg[p]) or pat_comment.search(arg[p]):
      arg.remove(arg[p])
      p -= 1
    else:
      p += 1
    n = len(arg)

  # Note: order matters

  # change `C code` to Ccode["C code"]
  n = len(arg)
  for i in range(n):
    if re.search(r"^\s*`", arg[i]) or re.search(r"`;?$", arg[i]):
      # check the syntax
      if not re.search(r"^\s*`", arg[i]) or not re.search(r"`;?", arg[i]):
        raise Exception(
            "The syntax of command {} is wrong, it must be between '`'.\n".
            format(arg[i]))
      arg[i] = re.sub(r"^\s*`", 'Ccode["', arg[i])
      arg[i] = re.sub(r"`;?", '"];', arg[i])

  # change "a(i) == b(i)" to Cpopulate["a=b"];
  for i in range(n):
    if (arg[i].find("==") != -1 and
        # exclude the following
        # NOTE: not all spaces removed yet!
        not re.search(r"(?i)file_name\s*=", arg[i]) and
        not re.search(r"(?i)Dimension\s*=", arg[i]) and
        not re.search(r"(?i)C_macro\d*\s*=", arg[i]) and
        not re.search(r"(?i)C_arg\d*\s*=", arg[i]) and
        not re.search(r"(?i)point\s*=", arg[i]) and
        not re.search(r"(?i)symm?\s*\[", arg[i]) and
        not re.search(r"(?i)comm?and?\s*\[", arg[i]) and
        not re.search(r"(?i)Pcode\s*\[", arg[i]) and
        not re.search(r"(?i)Declare\s*=", arg[i]) and
        not re.search(r"(?i)Ccode\s*\[", arg[i]) and
        not re.search(r"(?i)Cpopulate\s*\[", arg[i])):
      # split to lhs and rhs
      side = arg[i].split("==")
      # remove ; at rhs if any
      side[1] = re.sub(r";", "", side[1])
      # remove \n at rhs if any
      side[1] = re.sub(r"\n$", "", side[1])
      # remove (.+) at lhs (if any)
      side[0] = re.sub(r"\(.+\)", "", side[0])
      # remove (.+) at rhs (if any)
      side[1] = re.sub(r"\(.+\)", "", side[1])
      arg[i] = "Cpopulate[" + side[0] + "=" + side[1] + "];"

  # remove white spaces and comment at the middle
  # substitute ';' at the end with '@'
  n = len(arg)
  for i in range(n):
    if re.search(r"(?i)Ccode", arg[i]):
      arg[i] = re.sub(r"\s+$", "", arg[i])
      try:
        Ccode_part = re.search(r"(?i)Ccode\[.*\];?", arg[i]).group(0)
      except:
        raise Exception("Something seems wrong in {}.".format(arg[i]))
      parts = arg[i].split(Ccode_part)
      parts[0] = re.sub(r"\s+", "", parts[0])
      parts[1] = re.sub(r"\s+", "", parts[1])
      parts[1] = re.sub(r"#+.*", "", parts[1])
      arg[i] = parts[0] + Ccode_part + parts[1]
      if not re.search(r";$", arg[i]):
        raise Exception("Command {} must be ended with ';'".format(arg[i]))
      else:
        arg[i] = re.sub(r";$", "@", arg[i])
      continue
    elif re.search(r"(?i)comm?and", arg[i]):
      arg[i] = re.sub(r"\s+$", "", arg[i])
      try:
        Command_part = re.search(r"(?i)comm?and\[.*\];?", arg[i]).group(0)
      except:
        raise Exception("Something seems wrong in {}.".format(arg[i]))
      parts = arg[i].split(Command_part)
      parts[0] = re.sub(r"\s+", "", parts[0])
      parts[1] = re.sub(r"\s+", "", parts[1])
      parts[1] = re.sub(r"#+.*", "", parts[1])
      arg[i] = parts[0] + Command_part + parts[1]
      if not re.search(r";$", arg[i]):
        raise Exception("Command {} must be ended with ';'".format(arg[i]))
      else:
        arg[i] = re.sub(r";$", "@", arg[i])
      continue
    elif re.search(r"(?i)^C_macro\d*", arg[i]):
      # remove white spaces only before the first occurrence of '='
      arg[i] = re.sub(r"(?i)^C_macro(\d*)\s*=", "C_macro\\1=", arg[i])
      # remove comments
      arg[i] = re.sub(r"(?:\w)*#+.*", "", arg[i])
      # remove trailing white spaces
      arg[i] = re.sub(r"\s+$", "", arg[i])
      # remove leading white spaces
      parts = arg[i].split("=", 1)
      parts[1] = re.sub(r"^\s+", "", parts[1])
      arg[i] = parts[0] + "=" + parts[1]
      # add ;
      if not re.search(r";$", arg[i]):
        raise Exception("Command {} must be ended with ';'".format(arg[i]))
      else:
        arg[i] = re.sub(r";$", "@", arg[i])
      continue
    elif re.search(r"(?i)Pcode", arg[i]):
      arg[i] = re.sub(r"\s+$", "", arg[i])
      Pcode_part = re.search(r"(?i)Pcode\[.*\];?", arg[i]).group(0)
      parts = arg[i].split(Pcode_part)
      parts[0] = re.sub(r"\s+", "", parts[0])
      parts[1] = re.sub(r"\s+", "", parts[1])
      parts[1] = re.sub(r"#+.*", "", parts[1])
      arg[i] = parts[0] + Pcode_part + parts[1]
      if not re.search(r";$", arg[i]):
        raise Exception("Command {} must be ended with ';'".format(arg[i]))
      else:
        arg[i] = re.sub(r";$", "@", arg[i])
      continue
    arg[i] = re.sub(r"(?:\w)*#+.*", "", arg[i])
    arg[i] = re.sub(r" +", "", arg[i])
    if re.search(r"\n$", arg[i]): # check for semicolon
      if not re.search(r"[=\}\{;]+\n$", arg[i]) and not re.search(
          r"^\n$", arg[i]):
        raise Exception("Semicolon is missing at {}.".format(arg[i]))
    arg[i] = re.sub(r"\s+", "", arg[i])
    arg[i] = re.sub(r";$", "@", arg[i])

  # join discontinuities in terms like name={...}
  db_l = []
  n = len(arg)
  j = 0
  for i in range(n - 1):
    if i <= j and j != 0:
      continue
    elif re.search(r"=$", arg[i]) and re.search(r"^\{", arg[i + 1]):
      par = arg[i]
      j = i + 1
      while not re.search(r"\}$", arg[j]) and j != n:
        par += arg[j]
        j += 1
      par += arg[j]
      db_l.append(par)
    else:
      db_l.append(arg[i])

  if j != n - 1:
    db_l.append(arg[n - 1])

  # replace ';' with '@' for the Declare environment if any remained
  n = len(db_l)
  for i in range(n):
    if re.search(r"(?i)Declare=\{.*\}", db_l[i]):
      db_l[i] = re.sub(r";", "@", db_l[i])

  if CPI__glob_pr_flg == 1:
    pr("-")
    print("Input file after trimming:")
    for i in db_l:
      print(i)
    pr("-")

  return db_l


# read the input file and clean it from comments and white spaces
def parse_inputs():
  global CPI__glob_pr_flg
  global CPI__glob_2dim_flg
  global CPI__glob_Cpi_version
  global CPI__glob_simplification_flag

  notes = '"Converting Mathematical Equations From sympy to C"'
  parser = argparse.ArgumentParser(description=notes)

  parser.add_argument(
      "Cpi_file",
      type=argparse.FileType("r"),
      help="The Cpi file that is to be converted to C code",
  )

  parser.add_argument(
      "--print",
      action="store",
      dest="print_flag",
      type=str,
      help="Activate the printing flag by Y(es) for debug purposes",
  )
  parser.add_argument(
      "--simplify",
      action="store",
      dest="simplify_flag",
      type=str,
      help=
      "Instruct the default simplification. The default is the built-in simplify function. "
      "Options = {off: no simplification.| "
      "factor: use factor().| "
      "mathematica: use Mathematica FullSimplify[].}.",
  )

  # parser.add_argument('--2dim'   , action = 'store',   dest="dim2_flag", type=str, help = 'skip 2-dimension manifold checks(Y/N)')

  parser.add_argument(
      "--version",
      action="version",
      version="%(prog)s {}".format(CPI__glob_Cpi_version),
  )

  args = parser.parse_args()

  # open file
  input = args.Cpi_file.readlines()

  try:
    file_name = re.search(r"[\w\.]+$", args.Cpi_file.name).group(0)
  except:
    raise Exception("file name {} could not be trimmed.".format(
        args.Cpi_file.name))

  input.insert(0, "file_name=" +
               file_name) # adding the name of the input file to the beginning

  if re.search(r"(?i)Y(es)?", str(args.print_flag)):
    CPI__glob_pr_flg = 1
  else:
    CPI__glob_pr_flg = 0

  if re.search(r"(?i)off", str(args.simplify_flag)):
    CPI__glob_simplification_flag = 0
  elif re.search(r"(?i)factor", str(args.simplify_flag)):
    CPI__glob_simplification_flag = 1
  elif re.search(r"(?i)mathematica", str(args.simplify_flag)):
    CPI__glob_simplification_flag = 2
  else: # defualt value
    CPI__glob_simplification_flag = -1

  # if re.search(r'(?i)Y(es)?',str(args.dim2_flag)):
  # CPI__glob_2dim_flg = 1
  # else:
  #  CPI__glob_2dim_flg = 0

  # close file
  args.Cpi_file.close()

  if CPI__glob_pr_flg:
    pr("-")
    print("Input file:")
    print(input)
    pr("-")

  ## remove white spaces and comments
  input = trim_inputs(input)

  ## some doctests
  checkup_input(input)

  return input


class Maths_Info:

  def __init__(self, arg):
    # note: the order is important
    self.dim_i = self.parse_dimension(arg)
    self.c_macro_d = self.parse_c_macro(arg)
    self.c_arg_d = self.parse_c_arg(arg)
    self.symmetry_ld = self.parse_symmetry(arg)
    self.symbols_ld = self.parse_symbols(arg)
    self.instruction_dd = self.parse_instructions(arg)
    self.complete_symbols()
    self.populate_components()
    self.point_symb = self.parse_point_symbol(arg)
    self.indices_s = self.parse_indices(arg)

  # getting the parameter that used to show points on manifold:
  def parse_point_symbol(self, arg):
    flg = 0
    v = 0
    for s in arg:
      if re.search(r"(?i)^points?=", s):
        val = s.split("=")
        val[1] = re.sub(r"@", "", val[1])
        v = val[1]
        flg = 1

    # if no point given and it is required
    if flg == 0:
      flg_fatal = 0
      for obj in self.symbols_ld:
        if obj["obj"] == "field":
          flg_fatal = 1
          break

      if flg_fatal == 1:
        raise Exception("Can't find the manifold point symbol.\n")

    return v

  # realizing all of the indices in tensors:
  def parse_indices(self, arg):
    indx_s = set()
    n = len(self.instruction_dd)
    for i in range(n):
      if self.instruction_dd[str(i)]["job"] == "calc":
        s = self.instruction_dd[str(i)]["calc"]
        for obj in self.symbols_ld:
          s_tmp = s
          # if it is an indexed obj in the expression and not sclare
          if ("rank" in obj.keys() and int(obj["rank"]) > 0 and
              re.search(r"\b{}\b".format(obj["name"]), s)):
            tensor = "\\b{}\\b\\([-\\w,]+\\)".format(obj["name"])
            piece = re.findall(r"{}".format(tensor), s_tmp)
            # for each match
            for p in piece:
              p = re.sub(r"{}\(".format(obj["name"]), "", p)
              p = re.sub(r"\)".format(obj["name"]), "", p)
              l = p.split(",")
              for i in l:
                # if there is minus sign
                if re.search(r"-", i):
                  i = re.sub(r"-", "", i)
                # check if there is no quantity with the same name as tensor indices
                for sym in self.symbols_ld:
                  if sym["name"] == i:
                    raise Exception(
                        "'{}' used as both tensor indices and variables.\n"
                        "Note: arguments of a functions can not have tensor indices"
                        .format(i))

                indx_s.add(i)
    # if we need more indices for higher rank tensors add more indices
    max = 0
    N = len(indx_s)
    indj = 0
    while max < N + 5:
      rep_flg = 1
      while rep_flg == 1:
        ex_index = "{}{}".format(CPI__glob_index_stem, indj)
        # check it does not exist
        rep_flg = 0
        for obj in self.symbols_ld:
          if obj["name"] == ex_index:
            rep_flg = 1
            break
        indj += 1

      max += 1
      indx_s.add(ex_index)

    return indx_s

  # realizing the other symbols in the equations
  def complete_symbols(self):
    n = len(self.instruction_dd)
    for i in range(n):
      if self.instruction_dd[str(i)]["job"] == "calc":
        realize_eq_py(self, self.instruction_dd[str(i)]["calc"])

    # also we might have a C function with some argument, let's add them
    sym = set()
    for i in range(n):
      if self.instruction_dd[str(i)]["job"] == "calc":
        for obj in self.symbols_ld:
          if obj["obj"] == "function" and re.search(
              r"\b{}\b".format(obj["name"]),
              self.instruction_dd[str(i)]["calc"],
          ):
            reg = "{}\\({{1}}[-\\w,]+\\)".format(obj["name"])
            if not re.search(
                r"{}".format(reg),
                self.instruction_dd[str(i)]["calc"],
            ):
              raise Exception(
                  "Could not find the arguments of \nfunction {} in \n{}"
                  "try to use simpler argumnet for the functions".format(
                      obj["name"],
                      self.instruction_dd[str(i)]["calc"],
                  ))

            args = re.findall(
                r"{}".format(reg),
                self.instruction_dd[str(i)]["calc"],
            )
            for arg in args:
              s = re.sub(r"^{}\(".format(obj["name"]), "", arg)
              s = re.sub(r"\)$", "", s)
              ind = s.split(",")
              for _ in ind:
                _ = re.sub(r"^-", "", _)
                sym.add(_)

    for s in sym:
      # check if the args are already defined
      found = 0
      for obj in self.symbols_ld:
        if obj["name"] == s:
          found = 1

      # add to the data base
      d = dict()
      if found == 0:
        d["name"] = s
        d["obj"] = "variable"
        d["Ccall"] = "none"
        self.symbols_ld.append(d)

    # check some reserved name to not confilict:
    reserved = [] # dir(sympy)
    reserved.append("L")
    reserved.append("Tensors_db")
    reserved.append("KD")
    reserved.append("KD0i")
    reserved.append("KD1i")
    reserved.append("KD2i")
    reserved.append("KD3i")
    reserved.append("EIJK")
    reserved.append("job")
    reserved.append("C_instructions")

    for obj in self.symbols_ld:
      if obj["name"] in reserved:
        raise Exception('Name "{}" is reserved, use another name.\n'.format(
            obj["name"]))

    # add Kronecker delta and Levi-Chevita tensor to the CPI__db

    # Kronecker delta:
    d = dict()
    d["name"] = "KD"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = "2"
    self.symbols_ld.append(d)

    # Levi Cevita:
    d = dict()
    d["name"] = "EIJK"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = str(self.dim_i)
    self.symbols_ld.append(d)

    # Kronecker delta KD(0,i):
    d = dict()
    d["name"] = "KD0i"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = "1"
    self.symbols_ld.append(d)

    # Kronecker delta KD(1,i):
    d = dict()
    d["name"] = "KD1i"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = "1"
    self.symbols_ld.append(d)

    # Kronecker delta KD(2,i):
    d = dict()
    d["name"] = "KD2i"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = "1"
    self.symbols_ld.append(d)

    # Kronecker delta KD(3,i):
    d = dict()
    d["name"] = "KD3i"
    d["obj"] = "local"
    d["Ccall"] = "none"
    d["type"] = "tensor"
    d["rank"] = "1"
    self.symbols_ld.append(d)

  # populating components of symbols
  def populate_components(self):
    # compound symmetry to symbols
    for obj in self.symbols_ld:
      symmtype = []
      obj["symmetry"] = []
      for symm in self.symmetry_ld:
        if symm["name"] == obj["name"]:
          symmtype.append(symm["symmetry"])

      # check for duplicates
      perms = set()
      for symm in symmtype:
        item = "{}".format(str(symm["perm"]))
        perms.add(item)
      if len(perms) != len(symmtype):
        raise Exception("Similar symmetries were found for '{}'.".format(
            obj["name"]))

      obj["symmetry"] = symmtype

    # populating
    for obj in self.symbols_ld:
      if obj["name"] == "KD":
        (
            obj["matrix_comp"],
            obj["array_comp"],
        ) = populate_Kronecker_Delta(self.dim_i)
      elif obj["name"] == "KD0i":
        (
            obj["matrix_comp"],
            obj["array_comp"],
        ) = populate_Kronecker_Delta_i(self.dim_i, 0)
      elif obj["name"] == "KD1i":
        (
            obj["matrix_comp"],
            obj["array_comp"],
        ) = populate_Kronecker_Delta_i(self.dim_i, 1)
      elif obj["name"] == "KD2i":
        (
            obj["matrix_comp"],
            obj["array_comp"],
        ) = populate_Kronecker_Delta_i(self.dim_i, 2)
      elif obj["name"] == "KD3i":
        (
            obj["matrix_comp"],
            obj["array_comp"],
        ) = populate_Kronecker_Delta_i(self.dim_i, 3)
      elif obj["name"] == "EIJK":
        obj["matrix_comp"], obj["array_comp"] = populate_Levi_Civita(self.dim_i)
      elif (obj["obj"] == "variable" or obj["obj"] == "field" or
            obj["obj"] == "local"):
        obj["matrix_comp"], obj["array_comp"] = populate_components(obj, self)
      elif obj["obj"] == "function":
        continue
      else:
        raise Exception("No Job!")

  # prints its contents:
  def pr(self):
    pr("-")
    print("dimension = {}\n".format(self.dim_i))
    pr("-")
    print("Macros:\n")
    print(self.c_macro_d)
    pr("-")
    print("C_args:\n")
    print(self.c_arg_d)
    pr("-")
    print("Symbols:\n")
    for i in self.symbols_ld:
      print(i)
    pr("-")
    print("Indices:\n")
    for i in self.indices_s:
      print(i)
    pr("-")
    print("symmetry:\n")
    for i in self.symmetry_ld:
      print(i)
    pr("-")
    print("Instructions:\n")
    for i in self.instruction_dd.items():
      print(i)
    pr("-")

  # definitions
  def parse_symmetry(self, arg):
    sym_ld = []
    for s in arg:
      if re.search(r"(?i)symm?\[", s):
        if not re.search(r"\]@$", s):
          raise Exception(
              "Symmetry command {} has not been written correctly.".format(s))
        symmetry = dict()
        s0 = re.sub(r"\]@$", "", s)
        s0 = re.sub(r"(?i)^symm?\[", "", s0)
        try:
          symmetry["name"] = re.search(r"^\w+", s0).group(0)
          s0_sides = s0.split("=")
          if not re.search(r"^[+-]?{}\(".format(symmetry["name"]), s0_sides[1]):
            raise Exception(
                "Symmetry command {} has not been written correctly.".format(s))
        except:
          raise Exception(
              "Symmetry command {} has not been written correctly.".format(s))

        symmetry["symmetry"] = realize_symmetry(s0)
        sym_ld.append(symmetry)

    return sym_ld

  # parsing c_macros:
  def parse_c_macro(self, arg):
    d = dict()
    v = 0
    for s in arg:
      if re.search(r"(?i)^C_macro\d*", s):
        val = s.split("=", 1)
        val[1] = re.sub(r"@", "", val[1])
        d[val[0]] = val[1]
    return d

  # parsing c_args:
  def parse_c_arg(self, arg):
    d = dict()
    v = 0
    for s in arg:
      if re.search(r"(?i)^C_arg\d*", s):
        s_arg = re.sub(r"(?i)^(C_arg\d*)(=)", "\\1`", s)
        val = s_arg.split("`")
        val[1] = re.sub(r"@", "", val[1])
        d[val[0]] = val[1]

    return d

  # parsing instructions
  def parse_instructions(self, arg):
    # remove file_name, Dimension,C_macros etc that are not related
    # to a job directly. what remains are Ccode and Calculation etc.
    todo = [] # todo list
    i = 0
    for s in arg:
      if (not re.search(r"(?i)^file_name=", s) and
          not re.search(r"(?i)^Dimension=", s) and
          not re.search(r"(?i)^C_macro\d*=", s) and
          not re.search(r"(?i)^C_arg\d*=", s) and
          not re.search(r"(?i)^point=", s) and
          not re.search(r"(?i)^symm?\[", s) and
          not re.search(r"(?i)^comm?and?\[", s) and not re.search(r":=", s)):
        todo.append(s)

    if not todo:
      raise Exception("There is no instruction in the inputfile!")

    instruct = dict()
    inst_n = 0 # instruction number
    # NOTE: the following declare_bin must match declare_bin at
    # parse_symbols func. this is an ad-hoc solution to have different
    # used declare.
    declare_bin = 0

    for s in todo:
      # take care of C codes
      if re.search(r"(?i)^Ccode\[", s):
        d = dict()
        Ccode = re.sub(r'(?i)^Ccode\["?', "", s)
        Ccode = re.sub(r'"?\]@$', "", Ccode)
        d["job"] = "Ccode"
        d["Ccode"] = Ccode
        instruct[str(inst_n)] = d
        inst_n += 1
      # take care of python codes
      elif re.search(r"(?i)^Pcode\[", s):
        d = dict()
        Pcode = re.sub(r'(?i)^Pcode\["?', "", s)
        Pcode = re.sub(r'"?\]@$', "", Pcode)
        d["job"] = "Pcode"
        d["Pcode"] = Pcode
        instruct[str(inst_n)] = d
        inst_n += 1
      # take care of populate commands
      elif re.search(r"(?i)^Cpopulate\[", s):
        d = dict()
        pop = re.sub(r"(?i)^Cpopulate\[", "", s)
        pop = re.sub(r"\]@$", "", pop)
        d["job"] = "Cpopulate"
        d["Cpopulate"] = pop
        instruct[str(inst_n)] = d
        inst_n += 1
      # take care of declaration
      elif re.search(r"(?i)^Declare=", s):
        d = dict()
        d["job"] = "Cdeclare"
        d["declare_bin"] = declare_bin
        declare_bin += 1
        instruct[str(inst_n)] = d
        inst_n += 1

      # take care of calculations
      else:
        d = dict()
        s = re.sub(r"@$", "", s)
        d = dict()
        d["job"] = "calc"
        d["calc"] = s

        if not "=" in s:
          raise Exception("It does not support expression {}.".format(s))

        instruct[str(inst_n)] = d
        inst_n += 1

    # make sure the user defined the heads and name of C functions:
    first = instruct["0"]
    if first["job"] != "Ccode":
      raise Exception(
          "\nThe input file must have at least a Ccode command(back-ticks).\n")

    return instruct

  # find dimension:
  def parse_dimension(self, arg):
    flg = 0
    v = 0
    for s in arg:
      if re.search(r"(?i)dimension", s):
        val = s.split("=")
        val[1] = re.sub(r"@", "", val[1])
        v = int(val[1])
        flg = 1

    if flg == 0:
      raise Exception("Can't find the dimention of the manifold")

    # I should remember why I put this restriction :)
    if v < 3 and CPI__glob_2dim_flg == 0:
      raise Exception("The dimension of manifold must be greater than 2.")

    return v

  # parsing index object, variables, fields etc/
  # for declared objects in Declare {}
  # its list of dictionaries each dictionary
  # hold info about the objects
  def parse_symbols(self, arg):
    db_list = []
    protect_name = [] # dir(sympy)

    # parse each object
    n = len(arg)
    declare_bin = 0
    for i in range(n):
      if re.search(r"(?i)Declare=", arg[i]):
        ps = re.sub(r"(?i)Declare=", "", arg[i])
        ps = re.sub(r"^\{", "", ps)
        ps = re.sub(r"\}$", "", ps)
        sub0 = ps.split("@")
        sub0.remove("")
        for sub1 in sub0:
          d = dict() # dictionary

          # trim
          if not re.search(r"^\(", sub1) or not re.search(r"\)$", sub1):
            raise Exception(
                "Each declaration must be written in parentheses.\n")
          sub1 = re.sub(r"^\(", "", sub1)
          sub1 = re.sub(r"\)$", "", sub1)

          # objects
          if re.search(r"(?i)obj", sub1):
            try:
              s = re.search(r"(?i)obj=\w+", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))
            mem = s.split("=")
            d["obj"] = mem[1].lower()

          # names
          if re.search(r"(?i)name", sub1):
            try: # syntax example: name = T0(i,-j,k,-l,-m) or name = S2.
              # NOTE: only - sign allowed.
              s = re.search(r"(?i)name=\w+(\([-\w,]+\))?", sub1).group(0)
            except:
              raise Exception("'{}' was not written correctly.".format(sub1))

            # check syntax that does not have both parentheses and rank
            if s.find("(") != -1 and re.search(r"(?i)rank", sub1):
              raise Exception("'{}' was not written correctly.".format(sub1))

            mem = s.split("=")
            # get the head, i.e., T0(i,-j,k,-l,-m) head = T0
            d["name"] = re.search(r"\w+", mem[1]).group(0)

            # pars tensor names that include rank, like:
            # name = T(i,-j,k,-l,-m) to name = T, rank = UDUDD
            if re.search(r"\(.+\)", mem[1]):
              # get parentheses contents:
              in_parenth = re.search(r"\(.+\)", mem[1]).group(0)
              in_parenth = re.sub(r"^\(", "", in_parenth)
              in_parenth = re.sub(r"\)$", "", in_parenth)
              # split commas
              all_commas = in_parenth.split(",")
              # realize each comma to the corresponding rank
              # - = D and no - = U.
              d["rank"] = []
              d["type"] = []
              counter = 0
              for c in all_commas:
                if c.find("-") != -1:
                  d["type"].append("D")
                else:
                  d["type"].append("U")
                counter += 1

              d["rank"] = "{}".format(counter)

            if d["name"] in protect_name:
              raise Exception(
                  "The name {} is protected for Sympy library, use another name."
                  .format(d["name"]))

          # ranks and types
          if re.search(r"(?i)rank", sub1):
            try:
              s = re.search(r"(?i)rank=\w+", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))

            d["rank"] = []
            d["type"] = []
            mem = s.split("=") # mem[0] = rank
            counter = 0
            type = mem[1]
            if re.search(r"(?i)(0|none|scalar)", type):
              d["type"] = "scalar"
              d["rank"] = "0"
            else:
              while re.search(r"(?i)[UD]+", type):
                if re.search(r"(?i)^U{1}?", type):
                  type = re.sub(r"(?i)^U{1}?", "", type)
                  d["type"].append("U")
                  counter += 1
                elif re.search(r"(?i)^D{1}?", type):
                  type = re.sub(r"(?i)^D{1}?", "", type)
                  d["type"].append("D")
                  counter += 1
                else:
                  raise Exception(
                      "Please define contravariant or invariant of the tensor {}!"
                      .format(d["name"]))

              d[mem[0]] = "{}".format(counter)

          # special C argument supersedes point
          if re.search(r"(?i)C_arg\d*", sub1):
            try:
              s = re.search(r"(?i)C_arg\d*", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))
            d["C_argument"] = s
            try:
              d[s] = self.c_arg_d[s]
            except:
              raise Exception("The is no defined {} in {}!".format(s, sub1))

          # how to call
          if re.search(r"(?i)C_macro\d*", sub1):
            try:
              s = re.search(r"(?i)C_macro\d*", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))
            d["Ccall"] = s
            try:
              d[s] = self.c_macro_d[s]
            except:
              raise Exception("The is no defined {} in {}!".format(s, sub1))

          if re.search(r"(?i)None", sub1):
            try:
              s = re.search(r"(?i)None", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))
            d["Ccall"] = s.lower()

          if re.search(r"(?i)Ccode", sub1):
            try:
              s = re.search(r"(?i)Ccode", sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))
            d["Ccall"] = s
            try:
              s = re.search(r'(?i)Ccode\[".*"\]', sub1).group(0)
            except:
              raise Exception("{} was not written correctly.".format(sub1))

            s = re.sub(r'(?i)^Ccode\["', "", s)
            s = re.sub(r'"\]$', "", s)
            d["Ccode"] = s

          assert d["obj"] in ["field", "variable", "function"]

          if d["obj"] == "field":
            if not "rank" in d.keys():
              d["rank"] = "0"
              d["type"] = "scalar"

          # add bin for call
          d["declare_bin"] = declare_bin

          db_list.append(d)
        # for each Declare add to bin
        declare_bin += 1

    return db_list


# Levi_Civita:
def populate_Levi_Civita(rank):
  tab = ""
  comp = "'{}"
  append = "'.format(Eijk("
  append2 = "["
  indx = 0
  code = ""

  for i in range(rank): # fors
    if i > 0:
      append += ","
      append2 += ","
      tab += "\t"
    append += "arg[{}]".format(i)
    append2 += "_{}".format(i)
    if i > 0:
      code += "{}comp{} = []\n".format(tab, i)
    code += "{}for _{} in range({}):\n".format(tab, i, rank)
    indx = i

  # the most inner loop
  tab += "\t"
  append += "))"
  append2 += "]"
  code += "{}arg = {}\n".format(tab, append2)
  code += "{}comp{} = {}{}\n".format(tab, indx + 1, comp, append)
  code += "{}array.append(comp{})\n".format(tab, indx + 1)
  code += "{}comp{}.append(comp{})\n".format(tab, indx, indx + 1)

  # going upward of the nested loop
  for i in range(rank - 1, 0, -1):
    tab = re.sub(r"\s{1}?$", "", tab)
    code += "{}comp{}.append(comp{})\n".format(tab, i - 1, i)

  comp0 = [] # the indexed_obj
  array = []
  exec(code)

  return comp0, array


# Kronecker delta:
def populate_Kronecker_Delta(N):
  array = []
  matrix = []

  for i in range(N):
    row = []
    for j in range(N):
      if i == j:
        row.append("1.")
        array.append("1.")
      else:
        array.append("0.")
        row.append("0.")

    matrix.append(row)

  return matrix, array


# Kronecker delta(?,i):
def populate_Kronecker_Delta_i(N, q):
  array = []
  matrix = []

  # if dim > 4
  if N > 4:
    raise Exception(
        "This function has not been developed for more than 4 dimension.\n"
        '      One can develop this by looking up the keyword "KD?i" and\n'
        "      add analogously for higher dimension.\n")

  row = []
  for i in range(N):
    if i == q:
      row.append("1.")
      array.append("1.")
    else:
      array.append("0.")
      row.append("0.")

  matrix = row
  # matrix.append(row)

  return matrix, array


# given list of indices and number of desired indices
# it returns a tuple of indices in string format
def indices_tuple(CPI__db, n):
  assert n <= len(CPI__db.indices_s)
  assert n != 0

  indices = []
  i = 0
  for ind in CPI__db.indices_s:
    if i == n:
      break
    if re.search(r"{}".format(CPI__glob_index_stem), ind):
      indices.append(ind)
      i += 1

  assert i == n

  tup = "("
  N = len(indices)
  for i in range(N - 1):
    tup += "{},".format(indices[i])
  tup += "{})".format(indices[N - 1])

  return tup


# a fancy print for introduction
def intro_print():
  print("""
           ..Y8888888888888              8888888888888888888888888
          .Y8                               oo               oo
         .Y8                                oo               oo
        .Y8                                 oo               oo
       .Y8                                  oo               oo
      .Y8                                   oo               oo
      Y8                                    oo               oo
      Y8                                    oo               oo
      .Y8                                   oo               oo
       .Y8                                  oo               oo
        .Y8                                 oo               oo
         .Y8                                oo               oo
          .Y8                               oo               oo
           ..Y8888888888888                 oo               oo """)
  print("\n")
  print("Convertor of Mathematical Equations to C programming Language")
  print("\n\n")


# simplifying the expressions CPI__expr.
# one can make it more complete and instruct it with flags or inputs
# from the input file. for now, it's quite rudimentary and simple.
# it seems that most of the time doit() function alone is satisfactory
# because it cancels or adds the terms which reduces the number of
# operations. now, one can call extra factor() to simplify even more,
# however, this might add unnecessary constant numbers
# into the expression.
def my_simplify(CPI__expr):
  if CPI__glob_simplification_flag == -1:
    return simplify(CPI__expr, ratio=1)

  elif CPI__glob_simplification_flag == 0:
    return CPI__expr

  elif CPI__glob_simplification_flag == 1:
    return factor(CPI__expr)

  elif CPI__glob_simplification_flag == 2:
    return simplify_mathematica(CPI__expr)

  else:
    raise Exception(
        "No such simplification {}!".format(CPI__glob_simplification_flag))


# simplify using Mathematica
def simplify_mathematica(cpi_expr):
  s_repl = "XXCPIXX" # weird srt to replace '_'

  cpi_expr = str(cpi_expr)

  # check no similar substr
  if cpi_expr.find(s_repl) != -1:
    raise Exception("{} must not be used for a variable name!".format(s_repl))

  # replace '_' with s_repl for Mathematica
  cpi_expr = cpi_expr.replace("_", s_repl)

  # make a temp file
  pid = os.getpid()
  mfile_name = ".mfile_temp_{}".format(pid)
  omfile_name = ".omfile_temp_{}".format(pid)
  m_expr = mathematica_code(cpi_expr)
  m_code = m_expr + ";\n"
  m_code += "FullSimplify[%,TimeConstraint->2000];\n"
  # m_code += "Simplify[%,TimeConstraint->2000];\n"
  #### m_code += "c = CForm[%];\n"
  m_code += "Print[%]" + "\n"

  # write into mathematica file
  mfile = open(mfile_name, "w")
  mfile.write(m_code)
  mfile.close()

  # run mathematica and output into another
  cmd = "math -run -noprompt < {} 1> {}".format(mfile_name, omfile_name)
  ret = os.system(cmd)

  # read math results.
  omfile = open(omfile_name, "r")
  m_expr = omfile.read()
  omfile.close()

  # delete the files
  cmd = "rm -rf {} {}".format(mfile_name, omfile_name)
  os.system(cmd)

  # replace '$' with '_' for Cpi
  m_expr = m_expr.replace(s_repl, "_")

  return m_expr


# replace built-in variables (cpi_name,cpi_head,cpi_index) in macro_line
# note the args must be str.
def C_macro_replace_built_in(macro_line, cpi_name, cpi_head, cpi_index, tab):
  # NOTE: macro_line gets replaced.
  # replace CPI_name with cpi_name.
  # older version does not have prefix 'CPI_'.
  # thus, if 'CPI_name' prefix exists first check this,
  # otherwise replace 'name' (if any).
  if macro_line.find("CPI_name") != -1:
    macro_line = re.sub(r"CPI_name", cpi_name, macro_line)
  elif macro_line.find("name") != -1:
    macro_line = re.sub(r"name", cpi_name, macro_line)

  # replace CPI_index with cpi_index value.
  macro_line = re.sub(r"CPI_index", cpi_index, macro_line)

  # replace CPI_head with head value.
  macro_line = re.sub(r"CPI_head", cpi_head, macro_line)

  macro_line = tab + macro_line + "\n"

  return macro_line


# style the equation prints when written in the output file.
# for instance, a long streq is broken into few lines with some max width.
def style_eq_pr(streq, tab):
  styled = re.sub(r"(.{50,70}[^eE][\+\-\*/]+)\s*", "\\1\n", streq)
  styled += "\n"

  return styled


# given a tensor components and its symmetry,
# it reduces the number of components according to the given symmetry
# info about name, rank, etc is in obj
# this produces a dynamic code, for instance, as a reference:
#
# for _0 in range(2):
# 	matrix1 = []
# 	for _1 in range(2):
# 		arg = [_0,_1]
# 		perm = [0,1]
# 		sign = "-"
# 		if(arg[perm[0]] == arg[perm[1]]):
# 			if(sign == "-"):
# 				matrix2 = "0."
# 			else:
# 				matrix2 = indexed_obj_db[arg[0]][arg[1]]
# 		elif(arg[perm[0]] > arg[perm[1]]):
# 			h            = arg[perm[1]]
# 			arg[perm[1]] = arg[perm[0]]
# 			arg[perm[0]] = h
# 			if(sign == "-"):
# 				if indexed_obj_db[arg[0]][arg[1]].find("-",0,1) != -1:
# 					matrix2 = re.sub(r"^-","",indexed_obj_db[arg[0]][arg[1]])
# 				else:
# 					matrix2 = "-"+indexed_obj_db[arg[0]][arg[1]]
# 			else:
# 				matrix2 = indexed_obj_db[arg[0]][arg[1]]
# 		else:
# 			matrix2 = indexed_obj_db[arg[0]][arg[1]]
# 		matrix2 = re.sub(r"^-0.","0.",matrix2)
# 		array0.append(matrix2)
# 		matrix1.append(matrix2)
# 	matrix0.append(matrix1)
#
def reduce_symm_components(indexed_obj_db, symm, obj, CPI__db):
  matrix0 = [] # -> shared between the exec and this function
  array0 = [] # -> shared between the exec and this function
  rank = int(obj["rank"])

  ## write dynamic code:
  code = ""
  tab = "" # count number of tab for the dynamic code
  append = "" # ==> ex: [arg[0]][arg[1]]
  append2 = "[" # ==> ex: arg = [_0,_1]
  indx = 0
  for i in range(rank): # fors
    if i > 0:
      append2 += ","
      tab += "\t"
    append += "[arg[{}]]".format(i)
    append2 += "_{}".format(i)
    if i > 0:
      code += "{}matrix{} = []\n".format(tab, i)
    code += "{}for _{} in range({}):\n".format(tab, i, CPI__db.dim_i)
    indx = i

  tab += "\t"
  append2 += "]"
  code += "{}arg = {}\n".format(tab, append2)
  perm = symm["perm"]
  sign = symm["sign"]
  code += "{}perm = [{},{}]\n".format(tab, perm[0], perm[1])
  code += '{}sign = "{}"\n'.format(tab, sign)
  code += "{}if(arg[perm[0]] == arg[perm[1]]):\n".format(tab)
  code += '{}\tif(sign == "-"):\n'.format(tab)
  # NOTE: this 0. becomes important later when using the entries
  # thus, 0. MUST be at the beginning of the str with EXACTLY same syntax.
  code += '{}\t\tmatrix{} = "0."\n'.format(tab, indx + 1)
  code += "{}\telse:\n".format(tab)
  code += "{}\t\tmatrix{} = indexed_obj_db{}\n".format(tab, indx + 1, append)
  code += "{}elif(arg[perm[0]] > arg[perm[1]]):\n".format(tab)
  code += "{}\th            = arg[perm[1]]\n".format(tab)
  code += "{}\targ[perm[1]] = arg[perm[0]]\n".format(tab)
  code += "{}\targ[perm[0]] = h\n".format(tab)
  code += '{}\tif(sign == "-"):\n'.format(tab)
  # if it has already a minus sign at the begging remove the - sign
  # NOTE: this minus sign becomes important later when using the entries
  # thus, the MINUS SIGN MUST be at the beginning of the str with EXACTLY same syntax.
  code += '{}\t\tif indexed_obj_db{}.find("-",0,1) != -1:\n'.format(tab, append)
  code += '{}\t\t\tmatrix{} = re.sub(r"^-","",indexed_obj_db{})\n'.format(
      tab, indx + 1, append)
  code += "{}\t\telse:\n".format(tab)
  code += '{}\t\t\tmatrix{} = "-"+indexed_obj_db{}\n'.format(
      tab, indx + 1, append)
  code += "{}\telse:\n".format(tab)
  code += "{}\t\tmatrix{} = indexed_obj_db{}\n".format(tab, indx + 1, append)
  code += "{}else:\n".format(tab)
  code += "{}\tmatrix{} = indexed_obj_db{}\n".format(tab, indx + 1, append)
  # remove - from '-0.'
  code += '{0}matrix{1} = re.sub(r"^-0.","0.",matrix{1})\n'.format(
      tab, indx + 1)
  code += "{}array0.append(matrix{})\n".format(tab, indx + 1)
  code += "{}matrix{}.append(matrix{})\n".format(tab, indx, indx + 1)
  # going upward of the nested loop
  for i in range(rank - 1, 0, -1):
    tab = re.sub(r"\s{1}?$", "", tab)
    code += "{}matrix{}.append(matrix{})\n".format(tab, i - 1, i)

  try:
    exec(code)
  except:
    raise Exception("Symmetries belong to '{}' are wrong.".format(obj["name"]))

  return matrix0, array0


if __name__ == "__main__":
  main()