Graph(const vec1<vec1<VertexType> >& _points_in, int domain)
  { 
    vec1<vec1<VertexType> > _points = compressGraph(_points_in);
    if(_points.size() > domain)
      throw GAPException("Graph too large");
    edges = _points;
    edges.resize(domain);
    
    for(int i : range1(_points.size()))
    {
        int i_size = _points[i].size();
        for(int j = 1; j <= i_size; ++j)
        {
            if(_points[i][j].target() <= 0 || _points[i][j].target() > domain) {
                throw GAPException("Graph contains out-of-bounds vertex: " + toString(_points[i][j].target()));
            }
            
            if(_points[i][j].colour() < 0 ) {
                throw GAPException(" Graph contains invalid edge colour: " + toString(_points[i][j].colour()));
            }
            VertexType edge(i, _points[i][j].colour());
            if(directed)
            {
              edge = edge.flipped();
            }

            edges[_points[i][j].target()].push_back(edge);
        }
    }
    for(int i : range1(edges.size()))
    {
        std::set<VertexType> pntset(edges[i].begin(), edges[i].end());
        edges[i] = vec1<VertexType>(pntset.begin(), pntset.end());
    }
  }

 GAPStabChainWrapper getscc(const vec1<int>& v)
 {
     D_ASSERT(!fixed_base);
     GAP_callFunction(FunObj_ChangeStabChain, stabChain, GAP_make(v));
     GAPStabChainWrapper scc(stabChain);
     if(v.empty())
         return scc;
     // Now we have to walk down the stabilizer chain, until we find the position where v ends
     int pos = 1;
     while(true)
     {
         int orbit = scc.getOrbitStart();
         while(pos <= v.size() && v[pos] != orbit)
             pos++;
         if(pos > v.size())
         {
             return scc;
         }
         if(!scc.hasNextLevel())
             return scc;
         scc = scc.getNextLevel();
         if(!scc.hasOrbit())
         {
             // Reached bottom, so exit early. getOrbitsPartition still handles this case.
             return scc;
         }
     }
 }

std::string print_1d(const vec1<T>& vec){
  std::stringstream ss;
  ss << "[ ";
  for (size_t i = 0; i < vec.size()-1; ++i){
    ss << vec[i] << ", ";
  }
  ss << vec.back() << "]";
  return ss.str();
}

inline vec1<int> invertVecAsPermutation(const vec1<int>& v)
{
    vec1<int> ret(v.size());
    for(int i : range1(v.size()))
    {
        ret[v[i]] = i;
    }
    return ret;
}

    const vec1<OrbitalGraph>* fillRBaseOrbitalsCache(const vec1<int>& fix)
    {
        vec1<OrbitalGraph> orbitals = scc.orbitals(fix, ps->domainSize());

        if(original_orbitals.size() < fix.size() + 1)
            original_orbitals.resize(fix.size() + 1);
    
        original_orbitals[fix.size() + 1] = std::move(orbitals);

         return &(original_orbitals[fix.size() + 1]);
    }

    SplitState fix_buildingRBase(const vec1<int>& fixed_values, bool useOrbits, bool useBlocks, bool useOrbitals, bool rootCall = false)
    {
        debug_out(3, "scpg", "last depth " << last_depth.get() << " new depth " << fixed_values.size());
        D_ASSERT(fixed_values.size() > last_depth.get());
        last_depth.set(fixed_values.size());

        if(useOrbits)
        {
            doCacheCheck(config.useOrbits, tracking_first_found_orbits,
                         [this](const vec1<int>& v){ return this->fillRBaseOrbitPartitionCache(v); },
                         fixed_values, "orbits");
        }

        if(useBlocks)
        { 
             doCacheCheck(config.useBlocks, tracking_first_found_blocks,
                         [this](const vec1<int>& v){ return this->fillRBaseBlocksCache(v); },
                         fixed_values, "blocks");
        }

        if(useOrbitals)
        { 
            doCacheCheck(config.useOrbitals, tracking_first_found_orbitals,
                         [this](const vec1<int>& v){ return this->fillRBaseOrbitalsCache(v); },
                         fixed_values, "orbitals");
        }

        SplitState ss(true);

        int fixed_size = fixed_values.size();

        if(useOrbits)
        {
            const vec1<int>* part = 0;
            if(tracking_first_found_orbits.get() >= 0)
                part = this->getRBaseOrbitPartition_cached(tracking_first_found_orbits.get());
            else
                part = this->getRBaseOrbitPartition_cached(fixed_size);
            debug_out(3, "scpg", "fix_rBase:orbits");
            if(!part->empty())
                ss = filterPartitionStackByFunction(ps, SquareBrackToFunction(part));
            if(ss.hasFailed())
                return ss;
        }

        if( ( StabChainConfig::doConsiderIfNonTrivial(config.useOrbitals)
            && fixed_size == tracking_first_found_orbitals.get() ) ||
            ( config.useOrbitals == StabChainConfig::always ) ||
            ( rootCall ) )
        { return signal_changed_generic(range1(ps->cellCount()), identityPermutation()); }

        return ss;
    }

 const vec1<OrbitalGraph>* getRBaseOrbitals_cached(int s)
 {
     if(s < original_orbitals.size())
         return &(original_orbitals[s+1]);
     else
         return 0;
 }

 const vec1<std::map<int,int> >* getRBaseBlocks_cached(int s)
 {
     if(s < original_blocks.size())
         return &(original_blocks[s+1]);
     else
         return 0;
 }

 const vec1<int>* getRBaseOrbitPartition_cached(int s)
 {
     if(s < original_partitions.size())
         return &(original_partitions[s+1]);
     else
         return 0;
 }

    const vec1<std::map<int,int> >* fillRBaseBlocksCache(const vec1<int>& fix)
    {
        vec1<vec1<vec1<int> > > blocks = scc.blocks(fix);
        vec1<std::map<int,int> > block_functions;
        for(int i : range1(blocks.size()))
        {
            block_functions.push_back(partitionToMap(blocks[i]));
        }

        if(original_blocks.size() < fix.size() + 1)
            original_blocks.resize(fix.size() + 1);

        original_blocks[fix.size() + 1] = std::move(block_functions);

        return &(original_blocks[fix.size() + 1]);
    }

inline vec1<int> partitionToList(const vec1<vec1<int> >& part, int size, MissingPoints mp)
{
    vec1<int> vec(size, 0);
    int covered = 1;
    for(int i : range1(part.size()))
    {
        for(int val : part[i])
        {
            D_ASSERT(vec[val] == 0);
            vec[val] = i;
            covered++;
        }
    }    
    // Maybe the permutation is missing some values of the end. These
    // should all be treated as defined by 'MissingPoints'
    if(mp == MissingPoints_Fixed)
    {
        for(int i : range1(vec.size()))
        {
            if(vec[i] == 0)
                vec[i] = vec.size() + 1 + i; // (is +1 required? It doesn't hurt)
        }
    }
    return vec;
}

 bool check_transversal(const vec1<Permutation>& trans, int start, int end)
 {
     for(int i : range1(trans.size()))
     {
         start = trans[i][start];
     }
     return start == end;
 }

 vec1<int> getMins() const
 {
     vec1<int> minimums;
     for(int i : range1(orbit_mins.size()))
     {
         if(isMinimum(i))
             minimums.push_back(i);
     }
     return minimums;
 }

    void initalize(const vec1<int>& order)
    {
        D_ASSERT(!fixed_base);
        fixed_base = true;
        unpacked_stabChain_depth.resize(order.size());
        GAP_callFunction(FunObj_ChangeStabChain, stabChain, GAP_make(order));

        debug_out(1, "SCC", "Setting up cache");
        debug_out(3, "SCC", "Order " << order);
        int order_pos = 1;

        GAPStabChainWrapper stabChainCpy(stabChain);
        do
        {
            StabChainLevel scl(stabChainCpy);

            while(order[order_pos] != scl.base_value)
            {
                debug_out(3, "SCC", "Skipping depth " << order_pos);
                order_pos++;
            }

            debug_out(3, "SCC", "Setting depth "<<order_pos<<" base point "<<scl.base_value);
            levels.push_back(scl);
            unpacked_stabChain_depth[order_pos] = levels.size();

            stabChainCpy = stabChainCpy.getNextLevel();
        }
        while(stabChainCpy.hasNextLevel());

#ifndef NO_DEBUG
        for(int i : range1(unpacked_stabChain_depth.size()))
        {
            if(unpacked_stabChain_depth[i] != 0)
            {
                D_ASSERT(levels[unpacked_stabChain_depth[i]].base_value == order[i]);
            }
        }
#endif
    }

      HDINLINE float_X operator()(const float_X N, const float_X omega, const vec1 observer_unit_vec) const
      {

	/* Form Factor for CIC charge distribution of N discrete electrons:
	 * | \mathcal{F} |^2 = N + (N*N - N) * sinc^2(n_x * L_x * \omega) * sinc^2(n_y * L_y * \omega) * sinc^2(n_z * L_z * \omega)
	 *
	 * with observation direction (unit vector) \vec{n} = (n_x, n_y, n_z)
	 * and with: N     = weighting
	 *           omega = frequency
	 *           L_d   = the size of the CIC-particle / cell in dimension d
	 *
	 * the Form Factor: sqrt( | \mathcal{F} |^2 ) will be returned
	 */

	return sqrt(N + (N*N - N) * util::square(
						 math::sinc( observer_unit_vec.x() * CELL_WIDTH/(SPEED_OF_LIGHT*2)  * omega) *
						 math::sinc( observer_unit_vec.y() * CELL_HEIGHT/(SPEED_OF_LIGHT*2) * omega) *
						 math::sinc( observer_unit_vec.z() * CELL_DEPTH/(SPEED_OF_LIGHT*2)  * omega)
						  )
		    );

      }

// Takes a graph with multi-coloured (and possibly multiple occurrences)
// of edges, and replaces it with one where there is at most one edge between
// each pair of vertices.
inline vec1<vec1<ColEdge> > compressGraph(const vec1<vec1<ColEdge> >& graph)
{
  std::map<std::multiset<int>, int> seen_maps;
  
  vec1<vec1<ColEdge> > output_graph(graph.size());
  for(int i = 1; i <= graph.size(); i++) {
    std::map<int, std::multiset<int> > edges;
    for(int j = 1; j <= graph[i].size(); ++j) {
      edges[graph[i][j].target()].insert(graph[i][j].colour());
    }

    for(auto & edge : edges)
    {
      if(seen_maps.count(edge.second) == 0) {
        int val = seen_maps.size() + 1;
        seen_maps[edge.second] = val;
      }
      output_graph[i].push_back(ColEdge(edge.first, seen_maps[edge.second]));
    }
  }
  return output_graph;
}

    const vec1<int>* fillRBaseOrbitPartitionCache(const vec1<int>& fix)
    {
        debug_out(3, "scpg", "Fixing: "<< fix);
        vec1<vec1<int> > oart = scc.orbits(fix, ps->domainSize());
        debug_out(3, "scpg", "Got orbit partition" << oart);
        // This might not be necessary, but it doesn't hurt!
        for(int i : range1(oart.size()))
            std::sort(oart[i].begin(), oart[i].end());
        std::sort(oart.begin(), oart.end());

        vec1<int> filter;
        if(oart.size() > 1)
            filter = partitionToList(oart, ps->domainSize(), MissingPoints_Fixed);

        debug_out(3, "scpg", "Filter partition: "<< filter);

        if(original_partitions.size() < fix.size() + 1)
            original_partitions.resize(fix.size() + 1);

        original_partitions[fix.size()+1] = std::move(filter);

         return &(original_partitions[fix.size() + 1]);
    }

    void doCacheCheck(StabChainConfig::sc_config_option configchoice,
                      Reverting<int>& nontrivialdepth,
                      Func cache_fill,
                      const vec1<int>& fixed_values, const char* name)
    {
        (void)name;

        // Note that we always start with 'fixed_values' is empty, even if the group
        // contains some fixed points to start with
        if(configchoice == StabChainConfig::alwaysbase)
        {
            if(fixed_values.size() == 0 && nontrivialdepth.get() < 0)
            {
                const auto& ptr = cache_fill(fixed_values);
                if(ptr->size() > 0)
                {
                    nontrivialdepth.set(fixed_values.size());
                }
            }
        }
        else if(configchoice == StabChainConfig::firstnontrivial)
        {
            if(nontrivialdepth.get() < 0)
            {
                const auto& ptr = cache_fill(fixed_values);
                if(ptr->size() > 0)
                {
                    nontrivialdepth.set(fixed_values.size());
                }
            }
        }
        else
        {
            cache_fill(fixed_values);
        }
    }

inline std::map<int, int> partitionToMap(const vec1<vec1<int> >& part)
{
    std::map<int, int> m;
    int covered = 1;
    for(int i : range1(part.size()))
    {
        for(int val : part[i])
        {
            D_ASSERT(m[val] == 0);
            m[val] = i;
            covered++;
        }
    }    
    
    return m;
}

void broydn2(vec1<double>& x, T& func, const double eta,
             const double TOLF = 1e-8, int clear = 0) {
  typedef vec1<double> vec;
  const int n = x.size();

  vec F = func(x);
  if ((std::abs(*std::max_element(F.begin(), F.end(), [](double a, double b) {
        return std::abs(a) < std::abs(b);
      }))) <= 0.01 * TOLF) {
    std::cout << "Given initial guess is already a minimum\n";
  return;
  }

  vec xold = x;
  vec Fold = F;
  for (int i = 0; i < n; ++i){
    x[i] = xold[i] + eta * F[i];
  }

  //std::cout << eta << "\n";
  //std::cout << F << "\n";
  //std::cout << xold << x << "\n";


  int i = 0;
  bool check = false;
  vec2<double> u;
  vec2<double> v;

  while (not check) {
    F = func(x);
    vec Fdiff(F.size());
    for (int j = 0; j < n; ++j){
      Fdiff[j] = F[j] - Fold[j];
    }

    vec ui = x;
    for (int j = 0; j < n; ++j) {
      ui[j] -= xold[j] + eta * Fdiff[j];
    }

    for (int j = 0; j < i; ++j) {
      double scalar = 0;
      for (int l = 0; l < n; ++l){
        scalar += v[j][l] * Fdiff[l];
      }
      for (int l = 0; l < n; ++l) {
        ui[l] -=  scalar * u[j][l];
      }
    }
    u.push_back(ui);
    vec vi = Fdiff;
    double scalar = 0;
    for (int j = 0; j < n; ++j){
      scalar += Fdiff[j] * Fdiff[j];
    }
    if (scalar == 0){
      std::cerr << "norm in broydn2 is 0.0\n";
      std::exit(1);
    }
    for (int j = 0; j < n; ++j){
      vi[j] = Fdiff[j] / scalar;
    }
    v.push_back(vi);

    xold = x;
    Fold = F;
    for (int j = 0; j < n; ++j){
      x[j] = xold[j] + eta * F[j];
    }

    for (int j = 0; j < i+1; ++j){
      double scalar = 0;
      for (int l = 0; l < n; ++l){
        scalar += v[j][l] * F[l];
      }
      for (int l = 0; l < n; ++l){
        x[l] -=  scalar * u[j][l];
      }
    }
    //std::cout << F << "\n";
    double test =
        std::abs(*std::max_element(F.begin(), F.end(), [](double a, double b) {
          return std::abs(a) < std::abs(b);
        }));
    std::cout << i << ": " << test << "\n";
    if (test < TOLF){
      check = true;
    }

    //std::cout << x << "\n";
    i++;
    if (i == clear){
      u.clear();
      v.clear();
      i = 0;
      clear = 0;
    }
  }

//.........这里部分代码省略.........