using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices;
using System.IO;
using Graphs;
using Matching;
using MatrixLibrary;

namespace Matching
{
    unsafe public class Glpk_Solver : IGraphMatching
    {

        #region Attributes
        /*******************************************************************************************/
        /****************************************Attributes******************************************/
        /*******************************************************************************************/
        
        
        double* lp;
        Graph g1;
        Graph g2;
        Dictionary<Node, Node> nodeMatching;

        public Dictionary<Node, Node> NodeMatching
        {
            get { return nodeMatching; }
            set { nodeMatching = value; }
        }



        const string glpkLibrary = "glpk_4_49_86";
        readonly int GLP_FX = 5;
        readonly int GLP_DB = 4;       
        readonly int GLP_BV = 3;
        readonly int GLP_UNDEF=1;  /* solution is undefined */
        readonly int GLP_FEAS=2;   /* solution is feasible */
        readonly int GLP_INFEAS =3; /* solution is infeasible */
        readonly int GLP_NOFEAS =4; /* no feasible solution exists */
        readonly int GLP_OPT=5;     /* solution is optimal */
        readonly int GLP_UNBND=6;   /* solution is unbounded */

        readonly int GLP_ON = 1; /* enable something */
        readonly int GLP_OFF = 0;  /* disable something */






        struct ConstraintMatrix
        {
            public  int[] ia;
            public  int[] ja;
            public double[] ar;
        }
        #endregion

        #region DLLImport
        /*******************************************************************************************/
        /****************************************DLLImport******************************************/
        /*******************************************************************************************/

        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double* glp_create_prob();
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_add_rows(double* lp, int rows);
        [DllImport(glpkLibrary,CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_add_cols(double* lp, int cols);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_col_bnds(double* lp, int col, int bound_type, double lower_bound, double upper_bound);
        [DllImport(glpkLibrary,CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_col_kind(double* lp, int col, int kind);
        [DllImport(glpkLibrary,CallingConvention = CallingConvention.Cdecl,  SetLastError = true)]
        static extern void glp_set_obj_dir(double* lp, int dir);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_load_matrix(double* lp, int elements, int[] ia, int[] ja, double[] ar);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_col_name(double* P, int j, string name);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_row_name(double *P, int i, string name);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_row_bnds(double *P, int i, int type, double lb, double ub);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_set_obj_coef(double* P, int j, double coef);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_create_index(double *P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_find_col(double *P, string name);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_simplex(double *P, void * parm);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double glp_get_obj_val(double* P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double glp_get_col_prim(double* P, int j);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl,SetLastError = true,CharSet=CharSet.Ansi)]
        //[return: MarshalAs(UnmanagedType.LPStr)]
        static extern IntPtr glp_get_col_name(double *P, int j);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_get_num_cols(double * P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_delete_prob(double* lp);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_write_prob(double *P, int flags, string fname);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_print_sol(double *P, string fname);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_get_col_kind(double* P, int j);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_get_status(double * P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_get_num_int(double* P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_get_num_bin(double* P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_intopt(double *P, ref glp_iocp parm);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double glp_mip_obj_val(double * P);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double glp_mip_row_val(double* P, int i);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern double glp_mip_col_val(double * P, int j);
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern int glp_mip_status(double* P);    

        [StructLayout(LayoutKind.Sequential)]
        public struct glp_iocp
        { 
            /* integer optimizer control parameters */
            public  int msg_lev;            /* message level (see glp_smcp) */
            public int br_tech;            /* branching technique: */

            public int bt_tech;            /* backtracking technique: */
            public  int GLP_BT_DFS; /* depth first search */
            public  int GLP_BT_BFS ; /* breadth first search */
            public  int GLP_BT_BLB;  /* best local bound */
            public  int GLP_BT_BPH ; /* best projection heuristic */
            public double tol_int;         /* mip.tol_int */
            public double tol_obj;         /* mip.tol_obj */
            public int tm_lim;             /* mip.tm_lim (milliseconds) */
            public int out_frq;            /* mip.out_frq (milliseconds) */
            public int out_dly;            /* mip.out_dly (milliseconds) */
            //void (*cb_func)(glp_tree *T, void *info);
            [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
            public delegate void cb_func(void* T, void* info);
	   
            public void* cb_info;          /* mip.cb_info */
            public int cb_size;            /* mip.cb_size */
            public int pp_tech;            /* preprocessing technique: */
            public int GLP_PP_NONE  ; /* disable preprocessing */
            public   int GLP_PP_ROOT  ; /* preprocessing only on root level */
            public  int GLP_PP_ALL   ; /* preprocessing on all levels */
            double mip_gap;         /* relative MIP gap tolerance */
            public int mir_cuts;           /* MIR cuts       (GLP_ON/GLP_OFF) */
            public int gmi_cuts;           /* Gomory's cuts  (GLP_ON/GLP_OFF) */
            public int cov_cuts;           /* cover cuts     (GLP_ON/GLP_OFF) */
            public int clq_cuts;           /* clique cuts    (GLP_ON/GLP_OFF) */
            public int presolve;           /* enable/disable using MIP presolver */
            public int binarize;           /* try to binarize integer variables */
            public int fp_heur;            /* feasibility pump heuristic */
            public int alien;              /* use alien solver */
            public double[] foo_bar;     /* (reserved) */
         }
        
        [DllImport(glpkLibrary, CallingConvention = CallingConvention.Cdecl, SetLastError = true)]
        static extern void glp_init_iocp(ref glp_iocp parm);
        
        #endregion


        public Glpk_Solver()
        {
        }
        /// <summary>
        ///  constructor
        ///  Construct the problem and then solve it
        /// </summary>
        /// <param name="MainWindow"></param>
        public Glpk_Solver(Graph G1, Graph G2)
        {
            this.g1 = G1;
            this.g2 = G2;
            Console.WriteLine("Use the solver GLPK to solve the problem of Isomorphisme\n\n");
          
            int rownum=0;
            int colnum=0;
                        
            int numNodes1 = g1.ListNodes.Count();
            int numNodes2 = g2.ListNodes.Count();
            int numEdges1 = g1.ListEdges.Count();
            int numEdges2 = g2.ListEdges.Count();

           
            //for storing the value of the constraint matrix
            List<int> iaList = new List<int>();
            List<int> jaList = new List<int>();
            List<double> arList = new List<double>();
            iaList.Add(0);
            jaList.Add(0);
            arList.Add(0);
            lp = glp_create_prob();

            //construct le objet funcion
            for (int i = 1; i <= numNodes1; i++)
            {
                for (int k = 1; k <= numNodes2; k++)
                {
                    glp_add_cols(lp, 1);
                    glp_set_col_kind(lp, ++colnum, GLP_BV);
                    glp_set_col_name(lp, colnum, ("x_" + g1.ListNodes[i - 1].Id + "," + g2.ListNodes[k - 1].Id));
                   
                    double dist = g1.ListNodes[i - 1].Label.dissimilarity(g2.ListNodes[k - 1].Label);
                    glp_set_obj_coef(lp, colnum, dist);
                   
                }
            }
            
            for (int ij = 1; ij <= numEdges1; ij++)
                for (int kl = 1; kl <= numEdges2; kl++)
                {
                    glp_add_cols(lp, 1);
                    glp_set_col_name(lp, ++colnum, ("y_" + g1.ListEdges[ij - 1].Id + "," + g2.ListEdges[kl - 1].Id));
                    glp_set_col_kind(lp, colnum, GLP_BV);
                    glp_set_obj_coef(lp, colnum, g1.ListEdges[ij - 1].Label.dissimilarity(g2.ListEdges[kl - 1].Label));
                  
                }

            glp_create_index(lp);


            //construct contraints matrix
            ConstraintMatrix CM = new ConstraintMatrix();

            //contraint: equation 3 ( Every vertex of g1 must be matched to a unique vertex of g2)
            for (int i = 1; i <= numNodes1; i++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation3_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_FX, 1, 1);
                for (int k = 1; k <= numNodes2; k++)
                {
                    iaList.Add(rownum);
                    jaList.Add((i-1)*numNodes2+k);
                    arList.Add(1);
                }
                    

            }


            //equation 4 ( Every edge of g1 must be matched to a unique edge of g2)
            for (int ij = 1; ij <= numEdges1; ij++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation4_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_FX, 1, 1);
                for (int kl = 1; kl <= numEdges2; kl++)
                {

                    iaList.Add(rownum);
                    jaList.Add(numNodes1 * numNodes2 + (ij-1)*numEdges2+kl);
                    arList.Add(1);
                }
            }

            //equation 5 (Every vertex of g2 must be matched to at most an edge of g1 )
            for (int k = 1; k <= numNodes2; k++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation5_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_DB, 0, 1);
                for (int i = 1; i <= numNodes1; i++)
                {
                    iaList.Add(rownum);
                    jaList.Add(numNodes2 * (i - 1) + k);
                    arList.Add(1);
                }
            }

            //equation 6 (If two vertices are matched together, 
            //an edge originating one of these two vertices must be matched with an edge originating the other vertex)
            for (int k = 1; k <= numNodes2; k++)
                for (int ij = 1; ij <= numEdges1; ij++)
                {
                    glp_add_rows(lp, 1);
                    rownum++;
                    glp_set_row_bnds(lp, rownum, GLP_FX, 0, 0);
                    glp_set_row_name(lp, rownum, "equation6_" + rownum);
                    

                    string source = g1.ListEdges[ij-1].NodeSource.Id;
                    string nameVar = "x_" + source +  "," + g2.ListNodes[k - 1].Id;
                    int colInd = glp_find_col(lp, nameVar);
                    if(colInd==0) 
                        throw new InvalidProgramException();
                   
                    iaList.Add(rownum);
                    jaList.Add(colInd);
                    arList.Add(1);

                    string name1 = g1.ListEdges[ij - 1].Id;
                    foreach (Edge e in g2.ListNodes[k - 1].ListEdgesOut)
                    {
                        string name2 = e.Id;
                        nameVar = "y_" + name1 + "," + name2;
                        colInd = glp_find_col(lp, nameVar);
                        if (colInd == 0)
                            throw new InvalidProgramException();
                        
                        iaList.Add(rownum);
                        jaList.Add(colInd);
                        arList.Add(-1);
                    }
                }

            //equation 7 (If two vertices are matched together, 
            //an edge targeting one of these two vertices must be matched with an edge targeting the other vertex)
            for (int l = 1; l <= numNodes2; l++)
                for (int ij = 1; ij <= numEdges1; ij++)
                {
                    glp_add_rows(lp, 1);
                    glp_set_row_name(lp, ++rownum, "equation7_" + rownum);
                    glp_set_row_bnds(lp, rownum, GLP_FX, 0, 0);

                    string target = g1.ListEdges[ij - 1].NodeTarget.Id;
                    string nameVar = "x_" + target +","+ g2.ListNodes[l - 1].Id;
                    int colInd = glp_find_col(lp, nameVar);
                    if (colInd == 0)
                        throw new InvalidProgramException();

                    iaList.Add(rownum);
                    jaList.Add(colInd);
                    arList.Add(1);

                    string name1 = g1.ListEdges[ij - 1].Id;
                    foreach (Edge e in g2.ListNodes[l - 1].ListEdgesIn)
                    {
                        string name2 = e.Id;
                        nameVar = "y_" + name1 +","+ name2;
                        colInd = glp_find_col(lp, nameVar);
                        if (colInd == 0)
                            throw new InvalidProgramException();

                        iaList.Add(rownum);
                        jaList.Add(colInd);
                        arList.Add(-1);
                    }

                }

           
            int numEle = arList.Count()-1; 
            CM.ia = iaList.ToArray();
            CM.ja = jaList.ToArray();
            CM.ar = arList.ToArray();
            glp_load_matrix(lp, numEle, CM.ia, CM.ja, CM.ar);
            this.problemSave("prob.txt");
            solve_problem(lp);

            
                     
        }
     

        /// <summary>
        ///  destructor
        /// </summary>
        ~Glpk_Solver()
        {
            if(lp!=null)
            glp_delete_prob(lp);
        }

        /// <summary>
        ///  call the solve method to solve the problem transformed as a parameter
        /// </summary>
        /// <param name="MainWindow"></param>
        /// <param name="lp"></param>
        private void solve_problem(  double* lp)
        {
            //int num_int = glp_get_num_int(lp);
            //int num_bin = glp_get_num_bin(lp);

            glp_iocp param = new glp_iocp();
            param.GLP_BT_DFS = 1;
            param.GLP_BT_BFS = 2;  /* breadth first search */
            param.GLP_BT_BLB = 3;  /* best local bound */
            param.GLP_BT_BPH = 4;
            param.GLP_PP_NONE = 0; /* disable preprocessing */
            param.GLP_PP_ROOT = 1; /* preprocessing only on root level */
            param.GLP_PP_ALL = 2;
            param.foo_bar = new double[29];
            glp_init_iocp(ref param);

            param.presolve = 1;
            param.binarize = 1;

            int fal = glp_simplex(lp, null);
            int fal2 = glp_intopt(lp, ref param);
            print_sol(lp);
            this.solutionSave("sol.txt");
        }


        /// <summary>
        ///  print the result to the richText of the mainWindow
        /// </summary>
        /// <param name="MainWindow"></param>
        /// <param name="lp"></param>
        private void print_sol( double* lp)
        {
            //int status = glp_get_status(lp);
            int status = glp_mip_status(lp);

            if (status == GLP_NOFEAS)
            {
              
                Console.WriteLine("no feasible solution exists\n");
            }
            else if ( status == GLP_UNDEF)
                Console.WriteLine("solution is undefined \n");
            else if ( status == GLP_INFEAS)
                Console.WriteLine("solution is infeasible\n ");
            else if ( status == GLP_UNBND)
                Console.WriteLine("solution is unbounded\n ");
            else if (status==GLP_OPT||status==GLP_FEAS)
            {
                if (this.g1.ListNodes.Count() == this.g2.ListNodes.Count() && this.g1.ListEdges.Count() == this.g2.ListEdges.Count())
                    Console.WriteLine(" is the isomorphisme of " + this.g2.Name + "\n");
                else
                     Console.WriteLine(" is the sous-isomorphisme of " + this.g2.Name + "\n");
    

                double ress;
                double minimum = glp_mip_obj_val(lp);
                this.cost = minimum;
                Console.WriteLine("the minimum value of the objective function: " + minimum.ToString() + "\n");
                Console.WriteLine("///////////\n");

           
                this.nodeMatching = new Dictionary<Node,Node>();
                string valName;
                for (int i = 1; i <= (glp_get_num_cols(lp)); i++)
                {
                    //ress = glp_get_col_prim(lp, i);
                    ress = glp_mip_col_val(lp, i);
                    valName = System.Runtime.InteropServices.Marshal.PtrToStringAnsi(glp_get_col_name(lp, i));

                    if (ress == 1)
                    {
                        if (valName.StartsWith("x_"))
                        {
                            string[] res = seperateValName(valName);
                            Console.WriteLine("Nodes: " + this.g1.Name + "#" + res[0] + " ----------> " + this.g2.Name + "#" + res[1] + "\n");
                            Node source = g1.findNode(res[0]);
                            Node target = g1.findNode(res[1]);
                            nodeMatching.Add(source,target);
                        }
                        else if (valName.StartsWith("y_"))
                        {
                            string[] res = seperateValName(valName);
                            Console.WriteLine("Edges: " + g1.Name + "#" + res[0] + " ----------> " + g2.Name + "#" + res[1] + "\n");
                        }
                        else throw new InvalidProgramException("wrong name of the valeurs.");

                    }

                }
              

            }
        }

        /// <summary>
        ///  small tool function. Analysing the name of a variable
        ///  for exemple: valName="x_node1,node3", the function return a table of 2 elements: "node1", "node3"
        /// </summary>
        /// <param name="valName"></param>
        private string[] seperateValName(string valName)
        {
            string[] result;
            valName = valName.Substring(2);
            result = valName.Split(',');
            return result;
        }

        /// <summary>
        ///  save the problem in GLPK format into a text file.
        /// </summary>
        /// <param name="fname"></param>
        public int problemSave(string fname)
        {
            return glp_write_prob(this.lp, 0, fname);
        }

        /// <summary>
        ///  save the solve result into a text file.
        /// </summary>
        /// <param name="fname"></param>
        public int solutionSave(string fname)
        {
            return glp_print_sol(this.lp,fname);
        }

        public double cost { get; set; }

        public Dictionary<Node, Node> getNodeMatching()
        {
            throw new NotImplementedException();
        }

        public double getCost()
        {
            return this.cost;
        }

        public void Compute(Graph g1, Graph g2)
        {
            cost = -1;
            this.nodeMatching=new Dictionary<Node, Node>();
            this.g1 = g1;
            this.g2 = g2;
            Console.WriteLine("Use the solver GLPK to solve the problem of Isomorphisme\n\n");

            int rownum = 0;
            int colnum = 0;

            int numNodes1 = g1.ListNodes.Count();
            int numNodes2 = g2.ListNodes.Count();
            int numEdges1 = g1.ListEdges.Count();
            int numEdges2 = g2.ListEdges.Count();


            //for storing the value of the constraint matrix
            List<int> iaList = new List<int>();
            List<int> jaList = new List<int>();
            List<double> arList = new List<double>();
            iaList.Add(0);
            jaList.Add(0);
            arList.Add(0);
            lp = glp_create_prob();

            //construct le objet funcion
            for (int i = 1; i <= numNodes1; i++)
            {
                for (int k = 1; k <= numNodes2; k++)
                {
                    glp_add_cols(lp, 1);
                    glp_set_col_kind(lp, ++colnum, GLP_BV);
                    glp_set_col_name(lp, colnum, ("x_" + g1.ListNodes[i - 1].Id + "," + g2.ListNodes[k - 1].Id));

                    double dist = g1.ListNodes[i - 1].Label.dissimilarity(g2.ListNodes[k - 1].Label);
                    glp_set_obj_coef(lp, colnum, dist);

                }
            }

            for (int ij = 1; ij <= numEdges1; ij++)
                for (int kl = 1; kl <= numEdges2; kl++)
                {
                    glp_add_cols(lp, 1);
                    glp_set_col_name(lp, ++colnum, ("y_" + g1.ListEdges[ij - 1].Id + "," + g2.ListEdges[kl - 1].Id));
                    glp_set_col_kind(lp, colnum, GLP_BV);
                    glp_set_obj_coef(lp, colnum, g1.ListEdges[ij - 1].Label.dissimilarity(g2.ListEdges[kl - 1].Label));

                }

            glp_create_index(lp);


            //construct contraints matrix
            ConstraintMatrix CM = new ConstraintMatrix();

            //contraint: equation 3 ( Every vertex of g1 must be matched to a unique vertex of g2)
            for (int i = 1; i <= numNodes1; i++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation3_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_FX, 1, 1);
                for (int k = 1; k <= numNodes2; k++)
                {
                    iaList.Add(rownum);
                    jaList.Add((i - 1) * numNodes2 + k);
                    arList.Add(1);
                }


            }


            //equation 4 ( Every edge of g1 must be matched to a unique edge of g2)
            for (int ij = 1; ij <= numEdges1; ij++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation4_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_FX, 1, 1);
                for (int kl = 1; kl <= numEdges2; kl++)
                {

                    iaList.Add(rownum);
                    jaList.Add(numNodes1 * numNodes2 + (ij - 1) * numEdges2 + kl);
                    arList.Add(1);
                }
            }

            //equation 5 (Every vertex of g2 must be matched to at most an edge of g1 )
            for (int k = 1; k <= numNodes2; k++)
            {
                glp_add_rows(lp, 1);
                glp_set_row_name(lp, ++rownum, "equation5_" + rownum);
                glp_set_row_bnds(lp, rownum, GLP_DB, 0, 1);
                for (int i = 1; i <= numNodes1; i++)
                {
                    iaList.Add(rownum);
                    jaList.Add(numNodes2 * (i - 1) + k);
                    arList.Add(1);
                }
            }

            //equation 6 (If two vertices are matched together, 
            //an edge originating one of these two vertices must be matched with an edge originating the other vertex)
            for (int k = 1; k <= numNodes2; k++)
                for (int ij = 1; ij <= numEdges1; ij++)
                {
                    glp_add_rows(lp, 1);
                    rownum++;
                    glp_set_row_bnds(lp, rownum, GLP_FX, 0, 0);
                    glp_set_row_name(lp, rownum, "equation6_" + rownum);


                    string source = g1.ListEdges[ij - 1].NodeSource.Id;
                    string nameVar = "x_" + source + "," + g2.ListNodes[k - 1].Id;
                    int colInd = glp_find_col(lp, nameVar);
                    if (colInd == 0)
                        throw new InvalidProgramException();

                    iaList.Add(rownum);
                    jaList.Add(colInd);
                    arList.Add(1);

                    string name1 = g1.ListEdges[ij - 1].Id;
                    foreach (Edge e in g2.ListNodes[k - 1].ListEdgesOut)
                    {
                        string name2 = e.Id;
                        nameVar = "y_" + name1 + "," + name2;
                        colInd = glp_find_col(lp, nameVar);
                        if (colInd == 0)
                            throw new InvalidProgramException();

                        iaList.Add(rownum);
                        jaList.Add(colInd);
                        arList.Add(-1);
                    }
                }

            //equation 7 (If two vertices are matched together, 
            //an edge targeting one of these two vertices must be matched with an edge targeting the other vertex)
            for (int l = 1; l <= numNodes2; l++)
                for (int ij = 1; ij <= numEdges1; ij++)
                {
                    glp_add_rows(lp, 1);
                    glp_set_row_name(lp, ++rownum, "equation7_" + rownum);
                    glp_set_row_bnds(lp, rownum, GLP_FX, 0, 0);

                    string target = g1.ListEdges[ij - 1].NodeTarget.Id;
                    string nameVar = "x_" + target + "," + g2.ListNodes[l - 1].Id;
                    int colInd = glp_find_col(lp, nameVar);
                    if (colInd == 0)
                        throw new InvalidProgramException();

                    iaList.Add(rownum);
                    jaList.Add(colInd);
                    arList.Add(1);

                    string name1 = g1.ListEdges[ij - 1].Id;
                    foreach (Edge e in g2.ListNodes[l - 1].ListEdgesIn)
                    {
                        string name2 = e.Id;
                        nameVar = "y_" + name1 + "," + name2;
                        colInd = glp_find_col(lp, nameVar);
                        if (colInd == 0)
                            throw new InvalidProgramException();

                        iaList.Add(rownum);
                        jaList.Add(colInd);
                        arList.Add(-1);
                    }

                }


            int numEle = arList.Count() - 1;
            CM.ia = iaList.ToArray();
            CM.ja = jaList.ToArray();
            CM.ar = arList.ToArray();
            glp_load_matrix(lp, numEle, CM.ia, CM.ja, CM.ar);
            this.problemSave("prob.txt");
            solve_problem(lp);
        }
    }
}
