void compile_ex(const lst& exprs, const lst& syms, FUNCP_CUBA& fp, const std::string filename)
{
	lst replacements;
	for (std::size_t count=0; count<syms.nops(); ++count) {
		std::ostringstream s;
		s << "a[" << count << "]";
		replacements.append(syms.op(count) == symbol(s.str()));
	}

	std::vector<ex> expr_with_cname;
	for (std::size_t count=0; count<exprs.nops(); ++count) {
		expr_with_cname.push_back(exprs.op(count).subs(replacements));
	}

	std::ofstream ofs;
	std::string unique_filename = filename;
	global_excompiler.create_src_file(unique_filename, ofs);

	ofs << "void compiled_ex(const int* an, const double a[], const int* fn, double f[])" << std::endl;
	ofs << "{" << std::endl;
	for (std::size_t count=0; count<exprs.nops(); ++count) {
		ofs << "f[" << count << "] = ";
		expr_with_cname[count].print(GiNaC::print_csrc_double(ofs));
		ofs << ";" << std::endl;
	}
	ofs << "}" << std::endl;

	ofs.close();

	global_excompiler.compile_src_file(unique_filename, filename.empty());
	// This is not standard compliant! ... no conversion between
	// pointer-to-functions and pointer-to-objects ...
	fp = (FUNCP_CUBA) global_excompiler.link_so_file(unique_filename+".so", filename.empty());
}

matrix getJacobian(const lst& vf, const lst& expr, const lst& params){
	int i, j, nv, np;
	
	nv = vf.nops();
	np = params.nops();
	
	matrix Jac = matrix(nv, np);
	
	for (i = 0; i < nv; ++i)
	{
		ex f;
		if (expr == 0)
			f = vf[i];
		else
			f= iterated_subs(vf[i],expr);
		
		for (j = 0; j < np; ++j){
			symbol v = ex_to<symbol>(params[j]);
            ex df = f.diff(v);
            if (df != 0)
				Jac(i,j) = df;
		}
	}
	return Jac;
}

//
// Replace each var in e with delay(var,lag).
//
ex delay_vars(const ex& e, const ex& lag, const lst& vars)
{
    ex g = e;
    for (lst::const_iterator i = vars.begin(); i != vars.end(); ++i) {
        g = g.subs(*i == delay(*i, lag));
    }
    return g;
}

exset lst2set(const lst & l) {
  exset s;
  lst::const_iterator it = l.begin();
  while (it != l.end()) {
    s.insert(*it);
    ++it;
  }
  return s;
}

static void generate_lag_assignment(ofstream& fout, string name, const lst& lag,
        symbol& indvar, ex& history_formula)
{
    fout << "    if (" << indvar << " < " << lag.op(3) << ") {\n";
    fout << "        " << lag.op(0) << " = " << history_formula.subs(indvar == indvar - lag.op(3)) << endl;
    fout << "    }\n";
    fout << "    else {\n";
    fout << "        " << lag.op(0) << " = lagvalue(" << indvar - lag.op(3) <<  ", " <<  lag.op(1) << ")\n";
    fout << "    }\n";
}

/** Find all occurrences of a pattern. The found matches are appended to
 *  the "found" list. If the expression itself matches the pattern, the
 *  children are not further examined. This function returns true when any
 *  matches were found. */
bool ex::find(const ex & pattern, lst & found) const
{
	if (match(pattern)) {
		found.append(*this);
		found.sort();
		found.unique();
		return true;
	}
	bool any_found = false;
	for (size_t i=0; i<nops(); i++)
		if (op(i).find(pattern, found))
			any_found = true;
	return any_found;
}

void CloseContext(void)
{
 
	if (initialized) {

		for (lst::iterator i= MyHandleList.begin(); i != MyHandleList.end(); i++) {
			ThreadData* Tmp = *i;
			Tmp->thread_running = false;
				
		}
		Sleep(100);
		zmq_ctx_destroy (context);
		initialized = false;
	}
}

bool zero_volume(const lst& pole_list,const lst& w_list)
{
  try{
  using namespace lemon;
  GlpkLp lp;
  exhashmap<LpBase::Col> col_map;
  for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
    {
      col_map[*wi] = lp.addCol();
    }
  for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
    {
      ex tmp_expr = *pit;
      Lp::Expr constr_expr;
      for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
	{
	  
	  ex wi_coeff = tmp_expr.coeff(*wi);
	  tmp_expr-=(*wi)*wi_coeff;
	  if(is_a<numeric>(wi_coeff))
	    constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
	  else throw std::logic_error(std::string("Non numeric coefficient in pole term. "));
	}
      if(is_a<numeric>(tmp_expr))
	lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
      else 
	{
	  cout<< tmp_expr<<endl;
	  throw std::logic_error(std::string("Lower bound is not a numeric"));
	}
    }
    
  double l_bound,u_bound;
  for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
    {
      lp.min();
      lp.obj(col_map[*wi]);
      lp.solve();
      if (lp.primalType() == Lp::OPTIMAL) 
	l_bound = lp.primal();
      else return true;
      lp.max();
      lp.obj(col_map[*wi]);
      lp.solve();
      if (lp.primalType() == Lp::OPTIMAL) 
	u_bound = lp.primal();
      else return true;
      cout<<"ub: "<<u_bound<<" lb: "<<l_bound<<endl;
      if((u_bound - l_bound) <= 0)
        return true;
    }
  return false;
  }catch(std::exception &p)
    {
      throw std::logic_error(std::string("In function \"zero_volume\":\n |___> ")+p.what());
    }
}

ex iterated_subs(ex f, lst e)
{
    ex g = f;
    int n = e.nops();
    for (int i = n-1; i >= 0; --i) {
        g = g.subs(e[i]);
    }
    return g;
}

/** Substitute objects in an expression (syntactic substitution) and return
 *  the result as a new expression. */
ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
{
	GINAC_ASSERT(ls.nops() == lr.nops());

	// Convert the lists to a map
	exmap m;
	for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
		m.insert(std::make_pair(*its, *itr));

		// Search for products and powers in the expressions to be substituted
		// (for an optimization in expairseq::subs())
		if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
			options |= subs_options::pattern_is_product;
	}
	if (!(options & subs_options::pattern_is_product))
		options |= subs_options::pattern_is_not_product;

	return bp->subs(m, options);
}

string p_to_deriv_string(lst vars, int p[])
{
    int n = vars.nops();
    ostringstream dname;
    for (int j = 0; j < n; ++j) {
        for (int k = 0; k < p[j]; ++k) {
            dname << "_" << vars[j];
        }
    }
    return dname.str();
}

//
// Merge the expressions in exprs and formulas into a single expression.
// They are "merged" using the equality relation.  That is, if
// exprs = {e1, e2, e3} and formulas = {f1, f2}, the result is
// e1 == (e2 == (e3 == (f1 == f2)))
// Using == as the operator is just for convenience.  Any binary operation
// that does not occur in the expressions could have been used.
//
ex to_nested_tuple(const lst& exprs, const lst& formulas)
{
    lst all;

    for (lst::const_iterator iter = exprs.begin(); iter != exprs.end(); ++iter) {
        all.append(*iter);
    }
    for (lst::const_iterator iter = formulas.begin(); iter != formulas.end(); ++iter) {
        all.append(*iter);
    }

    int n = all.nops();
    assert(n != 0);

    if (n == 1) {
        return all.op(0);
    }
    ex t = all.op(n-2) == all.op(n-1);
    for (int k = n - 3; k >= 0; --k) {
        t = all.op(k) == t;
    }
    return t;
}

//
//  On return, lags is a GiNac::lst, where each element is a GiNac::lst
//  of length 4 containing {lagsym, variable_index + 1, var, lag_time}
//
void VectorField::convert_delay_to_lagvalue(ex& f, lst &lags)
{
    symbol t(IndependentVariable);
    exset dlist;
    f.find(delay(wild(1),wild(2)),dlist);
    for (exset::const_iterator iter = dlist.begin(); iter != dlist.end(); ++iter) {
        ex delayfunc = *iter;
        ex delayexpr = delayfunc.op(0);
        lst vars = FindVarsInEx(delayexpr);
        ex del = delayfunc.op(1);
        for (lst::const_iterator iter = vars.begin(); iter != vars.end(); ++iter) {
            ostringstream os;
            lst tmp;
            os << lags.nops() + 1;
            symbol lagsym("lag" + os.str());
            int vindex = FindVar(ex_to<symbol>(*iter));
            delayexpr = delayexpr.subs(*iter == lagsym);
            tmp = {lagsym, vindex + 1, *iter, del};
            lags.append(tmp);
        }
        f = f.subs(delayfunc == delayexpr);
    }
}

 GinacConstFunction::GinacConstFunction( const exvector& exs, const lst& params) 
     : Function( exs.size(), 0, params.nops()),
     exs_(exs), params_(params) {
         if(print__all){
         //std::cerr << "GinacConstFunction = ";					/**/
         //for(unsigned int i = 0; i < exs_.size(); i++ ){     /**/
         //    std::cerr << "[f("<<i<<") ="<<exs_[i] << "] " ; /**/
         //}													/**/
         //std::cerr << std::endl;								/**/
     	//}
     }
     
 void GinacConstFunction::eval(bool do_f, bool do_g, bool do_H){
     
     lst argsAndParams;
     lst::const_iterator params_it = params_.begin();
     
     for( unsigned int i = 0; params_it != params_.end(); params_it++, i++ ){
         argsAndParams.append( ex_to<symbol>(*params_it) == p(i));
     }
     //std::cerr << "NEXT FUNC"<< std::endl;           /**/
     
     if (do_f){
         for(unsigned int i = 0; i < exs_.size(); i++ ){
             f(i) = ex_to<numeric>(exs_[i].subs( argsAndParams )).to_double();
             //std::cerr << "f("<<i<<") ="<<f(i)<< std::endl; /**/
         }
     }
     
     if (do_g){
         clear_g();
     }
     
     if (do_H){
         clear_h();
     }
 }

UFXmap UF(lst k_lst,lst p_lst,lst subs_lst, unsigned int displacement )
{
  /// constructing expr Xi*Pi
  ex fprop = 0;
  ex sDtmp;
  string str;
  lst xp;
  for(lst::const_iterator it = p_lst.begin();it!=p_lst.end();++it)
    {
      str = "x_" + boost::lexical_cast<string>(displacement+std::distance(p_lst.begin(),it));
      xp.append(get_symbol(str)); // storing list of X identities
      fprop += get_symbol(str) * (*it); 
    }
  fprop = fprop.expand();
  sDtmp = fprop;
  //cout<<fprop<<endl;
  matrix M(k_lst.nops(),k_lst.nops());
  matrix Q(k_lst.nops(),1);//column
  GiNaC::numeric half(1,2);
  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    for(lst::const_iterator itc = itr;itc!=k_lst.end();++itc)
      if(itr == itc)
        {
          M(distance(k_lst.begin(),itr),distance(k_lst.begin(),itc)) = -1*fprop.coeff((*itr)*(*itc),1);
          //cout<<(*itr)*(*itc)<<"M("<<distance(k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr)*(*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr)*(*itc),1);
        }
      else
        {
          M(distance(k_lst.begin(),itr),distance( k_lst.begin(),itc)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          M(distance(k_lst.begin(),itc),distance(k_lst.begin(),itr)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          //cout<<(*itr)*(*itc)<<"M("<<distance( k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr),1).coeff((*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr),1).coeff((*itc),1);
        }
  cout<<"M: "<<M<<endl;
  sDtmp = sDtmp.expand();
  //cout<<"Expr linear on external momentum: "<<sDtmp.expand()<<endl;

  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    {
      Q(distance( k_lst.begin(),itr),0) = half*sDtmp.coeff((*itr),1);
      sDtmp -= (*itr)*sDtmp.coeff((*itr),1);
    }
  //  cout<<"Q: "<<Q<<endl;
  sDtmp = sDtmp.expand();
  ex minusJ = sDtmp;
   cout<<"-J: "<<minusJ<<endl;
  ex U = M.determinant();
  ex F = expand(M.determinant()*(minusJ+Q.transpose().mul(M.inverse()).mul(Q)(0,0)));
  lst lp;
  F=F.normal();
  cout<<"U= "<<U<<endl<<"F= "<<F<<endl;
  //  cout<<"pol_list "<<lp<<endl;
  U=U.subs(subs_lst,subs_options::algebraic);
  F=F.subs(subs_lst,subs_options::algebraic);
  cout<<"ALGEBRAIC: U= "<<U<<endl<<"F= "<<F<<endl;
  //  cout<<"UF_WORK"<<endl;
  return fusion::make_map<UFX::U,UFX::F,UFX::xlst>(U,F,xp);
}

/**
 *
 *  loop momentums,propagator expressions,
 *  invariants substitutions,propagator powers,number of loops
 *
 \param k_lst loop momentums list
 \param p_lst propagator expressions list
 \param subs_lst invariants substitutions list
 \param nu propagator powers list
 \param l number of loops 
   
 \return 
 *
 */
RoMB_loop_by_loop:: RoMB_loop_by_loop(
				      lst k_lst,
				      lst p_lst,
				      lst subs_lst,
                                      lst nu,
                                      bool subs_U
                                      )
{
  try
    {
      /* 
	 empty integral
      */
      MBintegral MBlbl_int(lst(),lst(),1);
      /* 
	 Full set of unused propagators, will change 
      */
      exlist input_prop_set;//( p_lst.begin(),p_lst.end());
      /* 
	 map for propagator powers
      */
      exmap prop_pow_map;
      for(lst::const_iterator Pit = p_lst.begin(); Pit != p_lst.end(); ++Pit)
	{
          input_prop_set.push_back(Pit->expand());
          prop_pow_map[Pit->expand()] = nu.op(std::distance(p_lst.begin(),Pit));
	}
	
      cout<<"INPSET: "<<input_prop_set<<endl;

      /* 
	 Iterate over momentums k1,k2,k3,etc.
      */
      unsigned int displacement_x = 0;
      unsigned int displacement_w = 0;
      for(lst::const_iterator kit = k_lst.begin(); kit != k_lst.end(); ++kit)
	{
          // Integral Normalization coefficient 
        //            MBlbl_int *= pow(I,k_lst.nops());
           MBlbl_int *= 1/tgamma(1+get_symbol("eps"));
          //MBlbl_int *= pow(Pi,2-get_symbol("eps"));
          //MBlbl_int *= exp(Euler*get_symbol("eps"));	  
          cout<<"PROP_POW_MAP "<<prop_pow_map<<endl;
	  /*
	    temporary set of propagators, with all momentum,except deleted
	  */
	  exlist tmp_p_lst(input_prop_set.begin(), input_prop_set.end()); 
	  /*
	    temporary set of propagators, with KIT momentum
	  */
	  lst P_with_k_lst;
	  BOOST_FOREACH(ex prop_tmp, tmp_p_lst)
	    {
	      if(prop_tmp.has(*kit))
		{
		  P_with_k_lst.append(prop_tmp);
		  input_prop_set.remove(prop_tmp);
		}
	    }
            
	  cout<< "Set wo k_i "<<input_prop_set<<endl;
	  cout<<" PWKlst "<<P_with_k_lst<<endl;
	  bool direct_formula_applied = false;
          // if only one term in PWKLST use well known formulas
          // [Smirnov A.1]
          if(!direct_formula_applied && (P_with_k_lst.nops() == 1))
            {
              ex pr_t = P_with_k_lst.op(0);
              ex nu_t = prop_pow_map[pr_t];
              exmap repls;
              BOOST_ASSERT_MSG(pr_t.match(-pow(*kit,2) + wild(2)),"ONE PROP");
              if(pr_t.match(-pow(*kit,2) + wild(2),repls))cout<<"repls: "<<repls<<endl;
              ex mass_tadpole = (tgamma(nu_t+get_symbol("eps")-2)/tgamma(nu_t)*pow(wild(2).subs(repls),-nu_t-get_symbol("eps")+2));
              cout<<mass_tadpole<<endl;
              MBlbl_int *= mass_tadpole;
              MBlbl_int.add_pole(nu_t+get_symbol("eps")-2);
              direct_formula_applied = true;
            }
          if(!direct_formula_applied && (P_with_k_lst.nops() == 2)) 
            {
              //TWO terms in PWK_LST, [Smirnov A.4]
              exmap repls_tad;
              if((P_with_k_lst.nops()==2) &&
        	 (( (P_with_k_lst.op(0).match(-pow(*kit,2))) && (P_with_k_lst.op(1).match(-pow(*kit,2)+wild())))||
                  ( (P_with_k_lst.op(1).match(-pow(*kit,2))) && (P_with_k_lst.op(0).match(-pow(*kit,2)+wild()))))
                 && !wild().has(*kit)
//.........这里部分代码省略.........

std::pair<ex,ex>  hyper_cube_den(lst pole_list,lst w_list, ex den)
{
try
  {
    using namespace lemon;
    GlpkLp lp;
    exhashmap<LpBase::Col> col_map;
    for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      {
        col_map[*wi] = lp.addCol();
      }
    for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
      {
        ex tmp_expr = *pit;
        // cout<<"Expr: "<<*pit<<" subexprs:"<<endl;
        Lp::Expr constr_expr;
        for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
          {
            // cout<<*wi<<"  coeff "<<tmp_expr.coeff(*wi)<<endl;
            ex wi_coeff = tmp_expr.coeff(*wi);
            tmp_expr-=(*wi)*wi_coeff;
            if(is_a<numeric>(wi_coeff))
              {
                constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
                //                cout<<ex_to<numeric>(wi_coeff).to_double()<<endl; 
              }
            else throw std::logic_error("Non numeric coefficient in pole term. ");
          }
        //constr_expr+=
        //cout<<"Ostatok "<<tmp_expr<<endl;
        if(is_a<numeric>(tmp_expr))
          lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
        else throw std::logic_error(std::string("Lower bound is not a numeric"));
      }

    double l_bound,r_bound;
    cout<<"Hyper cube"<<endl;
    //   for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      
        lp.min();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<*w_list.begin()<<" = ("<<lp.primal()<<",";
            l_bound = lp.primal();
          }
        else throw std::logic_error("Optimal solution not found.");
        lp.max();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<lp.primal()<<")"<<endl;
            r_bound = lp.primal();
          }
        else throw  std::logic_error("Optimal solution not found.");

        cout<<"Simplex : "<<endl;
        ex lbs,ubs;
        //        bool founds;
        //        tie(lbs,ubs,founds) = simplex_ara(pole_list,w_list,*w_list.begin());
        // cout<<"lbs: "<<lbs<<" ubs: "<<ubs<<endl;
        //  boost::mt19937 rng;       
        //        boost::uniform_real<> bounded_distribution(l_bound,r_bound);      // distribution that maps to 1..6
        // see random number distributions
        //  boost::variate_generator<boost::mt19937&, boost::uniform_real<> >
        //  die(rng, bounded_distribution);             // glues randomness with mapping
        //rng.seed(static_cast<unsigned char>(std::time(0)));

        
        //        double x   = die();                      // simulate rolling a die
        // new edition in center of interval no random
        //if((l_bound - trunc(l_bound)) < 10e-8) lbs = trunc(l_bound);
        //lbs = ginac_set_d(l_bound);
	lbs = l_bound;
        //if((r_bound - trunc(r_bound)) < 10e-8) ubs = trunc(r_bound);
	//ubs = ginac_set_d(r_bound);
	ubs = r_bound;
        //cout<<"ginac_set_d "<<ginac_set_d(l_bound)<<" "<<ginac_set_d(r_bound)<<endl;
        ex x_half = lbs*(1- pow(den,-1)) + ubs*pow(den,-1);
        //        double real_half = (r_bound + l_bound)/2.0;
        //        cout<<"Point "<<die()<<"   "<<x<<endl;
        return std::make_pair(*w_list.begin(),x_half);
  }catch(std::exception &p)
    {
      throw std::logic_error(std::string("In function \"hyper_cube_den\":\n |___> ")+p.what());
    }
}

void generate_deriv(string lang, ofstream &fout, ofstream &pout, string name, int r,
                         lst vf, lst expreqn_list,
                         lst vars, lst params)
{
    if (r < 1) {
        cerr << "Error: generate_deriv called with order=" << r << ". order must be greater than 0.\n";
        exit(-1);
    }
    //
    // Print the function header
    //
    fout << endl;
    fout << "/*\n";
    if (r == 1) {
        fout << " *     deriv = Df(x)[v1]\n";
    }
    else if (r == 2) {
        fout << " *     deriv = D^2f(x)[v1,v2]\n";
    }
    else if (r == 3) {
        fout << " *     deriv = D^3f(x)[v1,v2,v3]\n";
    }
    else {
        fout << " *     deriv = D^" << r << "f(x)[v1,...,v" << r << "]\n";
    }
    fout << " */\n";
    if (lang == "c") {
        fout << "void " << name << "_diff" << r << "(double deriv[], double x_[], double p_[],";
        pout << "void " << name << "_diff" << r << "(double deriv[], double x_[], double p_[],";
    }
    else {
        fout << "function " << name << "_diff" << r << "(x_, p_,";
    }
    for (int k = 0; k < r; ++k) {
        if (lang == "c") {
            fout << "double v" << k+1 << "_[]";
            pout << "double v" << k+1 << "_[]";
        }
        else {
            fout << "v" << k+1 << "_";
        }
        if (k < r-1) {
            fout << ",";
            if (lang == "c") {
                pout << ",";
            }
        }
    }
    fout << ")\n";
    if (lang == "c") {
        pout << ");\n";
    }
    fout << "{\n";
    if (lang == "c") {
        CDeclare(fout,"double",vars);
        CDeclare(fout,"double",params);
        fout << endl;
    }
    const char* pre = "    ";
    if (lang == "javascript") {
        pre = "    var ";
    }
    GetFromVector(fout, pre, vars, "=", "x_", "[]", 0, ";");
    fout << endl;
    GetFromVector(fout,pre,params, "=", "p_", "[]", 0, ";");
    fout << endl;
/*
    int na = exprnames.nops();
    for (int i = 0; i < na; ++i) {
        fout << "    " << exprnames[i] << " = " << exprformulas[i] << ";" << endl;
    }
    if (na > 0) {
        fout << endl;
    }
*/

    if (lang == "javascript") {
        fout << "    var deriv = [];\n";
    }
    int n = vars.nops();
    int *p = new int[n];
    for (int i = 0; i < n; ++i) {
        if (lang == "c") {
            fout << "{\n";
        }
        // ex f = vf[i];
        ex f = iterated_subs(vf[i],expreqn_list);
        map<string,ex> derivatives;
        // Initialize p to [r,0,0,...,0]
        p[0] = r;
        for (int k = 1; k < n; ++k) {
            p[k] = 0;
        }
        fout << "    /*\n";
        fout << "     *  Partial derivatives of vf[" << i << "]. \n";
        fout << "     *  Any derivative not listed here is zero.\n";
        fout << "     */\n";
        do {
            // Compute the partial derivative
            ex df = f;
//.........这里部分代码省略.........

void mexFunction( int nlhs, mxArray *plhs[], 
				 int nrhs, const mxArray *prhs[] ) 
{

   if(!initialized) initialize_zmq();
    
	if (nrhs < 2) {
		return;
	}


	char *Command = mxArrayToString(prhs[0]);

	if   (strcmp(Command, "StartConnectThread") == 0)  {
        	int url_length = int(mxGetNumberOfElements(prhs[1])) + 1;
			char *url = mxArrayToString(prhs[1]);
            ThreadData* sd=  create_client_thread(url, url_length);
			MyHandleList.push_back(sd);
	        plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL);
			double* Tmp= (double*)mxGetPr(plhs[0]);
			memcpy(Tmp, &sd, 8);
    }
    
    
    if   (strcmp(Command, "Send") == 0)  {
		double* Tmp= (double*)mxGetPr(prhs[1]);
		ThreadData* sd;
        memcpy(&sd, Tmp, 8);

		if (!sd->thread_running) {
			// Thread was killed and now we try to send things again?
			mexPrintf("Connection was closed. Cannot send.\n");
			return;
		}

		char *s= mxArrayToString(prhs[2]);
		std::string std_s(s);
		MyQueue.push(std_s);

	}
			
    if   (strcmp(Command, "CloseThread") == 0)  {
		double* Tmp= (double*)mxGetPr(prhs[1]);
		ThreadData* sd;
        memcpy(&sd, Tmp, 8);
		sd->thread_running = false;

		// remove from active handle list
		
		MyHandleList.remove(sd);
			
		//MyHandleList.push_back(sd);

	}


}