This is a wonderful Topcoder problem that draws from material that we've learned about Depth-First-Search and enumeration.

The first thing to do is to wrap your head around the problem through the examples. There are two graphs here, both with cities as nodes. The first is a complete graph, because there is a road from each city to each other city. The second is composed only of the edges representing roads that you have to take. These are specified by the adjacency matrix roads.

Let's go through the examples. For each example, I'll have the adjacency matrix on the left, and the graphs on the right. The graph from roads will have red edges. The other edges in the complete graph will be in gray.

Example 0:

{"NYN",
 "YNN",
 "NNN"}

Answer = 4

Example 1:

{"NYYY",
 "YNNN",
 "YNNN",
 "YNNN"}

Answer = 0

Example 2:

{"NYY",
 "YNY",
 "YYN"}
Answer = 0

Example 3:

{"NNNNNY",
 "NNNNYN",
 "NNNNYN",
 "NNNNNN",
 "NYYNNN",
 "YNNNNN"}

Answer = 24

Since you must visit all N cities and you must use exactly N-1 roads, you must take a simple path that visits all cities. For that reason Examples 1 and 2 return zero, because there's no way to construct a simple path with those graphs. Put another way, when you consider the red edges, every node must have two or fewer edges on its adjacency list, and there can be no cycles. So, our first programming task is going to be to create the graph with red edges, and to return 0 when we discover cases like examples 1 and 2.

That code is in Hamilton-1.cpp. The cycle detection should look very familiar to you from the lecture notes on Depth-First-Search: 

#include <string>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstdlib>
using namespace std;

class Vertex {
  public:
    int visited;
    vector <int> adj;
};

class HamiltonPath {
  public:
    int Cycle_Det(int n, int from);
    vector <Vertex> V;
    int countPaths(vector <string> roads);
};

int HamiltonPath::Cycle_Det(int n, int from)
{
  int i, j;

  if (V[n].visited != -1) return 1;

  V[n].visited = 1;
  for (i = 0; i < V[n].adj.size(); i++) {
    j = V[n].adj[i];
    if (j != from && Cycle_Det(j, n)) return 1;
  }
  return 0;
}

int HamiltonPath::countPaths(vector <string> roads)
{
  int i, j;

  V.resize(roads.size());
  for (i = 0; i < V.size(); i++) {
    V[i].visited = -1;
  }

  for (i = 0; i < roads.size(); i++) {     // Create the graph
    for (j = 0; j < roads.size(); j++) {
      if (roads[i][j] == 'Y') V[i].adj.push_back(j);
    }
  }

  /* Return 0 if a node has more than two edges coming from it. */
  for (i = 0; i < V.size(); i++) if (V[i].adj.size() > 2) return 0;

  /* Run Cycle Detection and return 0 if you discover a cycle */

  for (i = 0; i < V.size(); i++) {
    if (V[i].visited == -1 && Cycle_Det(i, -1)) return 0;
  }

  /* Return 1 for now -- we'll finish this later. */
  return 1;
}

To compile this, the file Hamilton-Main.cpp includes HamiltonPath.cpp and runs the examples from the command line. So we copy Hamilton-1.cpp to HamiltonPath.cpp and then compile Hamilton-Main.cpp. Examples 1 and 2 return 0 as expected, while the other two examples return 1: 

UNIX> cp Hamilton-1.cpp HamiltonPath.cpp
UNIX> g++ Hamilton-Main.cpp
UNIX> a.out 0
1
UNIX> a.out 1
0
UNIX> a.out 2
0
UNIX> a.out 3
1
UNIX>
At this point, we have to figure out how to count the paths. Let's consider examples. I'm only drawing the red edges now. Here's a pretty simple example:

A path through the graph is equivalent to a permutation of the numbers 0, 1, 2, 3. The number of paths is therefore (4!) = 24. How about example 0 from the writeup:

Well, there are two connected components, and every path through the graph will correspond to a permutation of the components. However, for each permutation of components, there are two ways to travel through nodes 0 and 1. Thus, the number of paths is (2!)2 = 4. Finally, let's look at example 3 from the writeup:

Again, each path corresponds to a permutation of the three components. For the first two components, there are two ways to go through the component. Therefore, the number of paths is (3!)(2)(2) = 24.

Do you see the pattern? If there are c connected components and of these c, there are b components that have two or more nodes, then the number of paths through the graph is:

(c!) 2b

So, to solve this problem, let's identify the connected components, which will give us c and b. Then we can calculate the above formula, taking the product modulo 1,000,000,007 at each step. It's best to use a 64-bit integer for these calculations. That's why you'll see the term long long. The solution is in Hamilton-2.cpp 

#include <string>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstdlib>
using namespace std;

class Vertex {
  public:
    int visited;
    int component;
    vector <int> adj;
};
class HamiltonPath {
  public:
    void ConComp(int n, int c);
    int Cycle_Det(int n, int from);
    vector <Vertex> V;
    int countPaths(vector <string> roads);
};

void HamiltonPath::ConComp(int n, int c)
{
  int i;

  if (V[n].component != -1) return;
  V[n].component = c;
  for (i = 0; i < V[n].adj.size(); i++) {
    ConComp(V[n].adj[i], c);
  }
}

int HamiltonPath::Cycle_Det(int n, int from)
{
  int i, j;

  if (V[n].visited != -1) return 1;

  V[n].visited = 1;
  for (i = 0; i < V[n].adj.size(); i++) {
    j = V[n].adj[i];
    if (j != from && Cycle_Det(j, n)) return 1;
  }
  return 0;
}

int HamiltonPath::countPaths(vector <string> roads)
{
  int i, j;
  int c, b;
  long long p;

  V.resize(roads.size());
  for (i = 0; i < V.size(); i++) {
    V[i].visited = -1;
    V[i].component = -1;
  }

  for (i = 0; i < roads.size(); i++) {
    for (j = 0; j < i; j++) {
      if (roads[i][j] == 'Y') {
        V[i].adj.push_back(j);
        V[j].adj.push_back(i);
      }
    }
  }
  for (i = 0; i < V.size(); i++) if (V[i].adj.size() > 2) return 0;
  for (i = 0; i < V.size(); i++) {
    if (V[i].visited == -1 && Cycle_Det(i, -1)) return 0;
  }

  /* Determine connected components */
  c = 0;
  b = 0;
  for (i = 0; i < V.size(); i++) {
    if (V[i].component == -1) {
      if (V[i].adj.size() > 0) b++;
      ConComp(i, c);
      c++;
    }
  }

  /* Calculate (c!)2^b */
  p = 1;
  for (i = 2; i <= c; i++) {
    p *= i;
    p %= 1000000007;
  }
  for (i = 0; i < b; i++) {
    p *= 2;
    p %= 1000000007;
  }
  return p;
}

Compile and test: 

UNIX> cp Hamilton-2.cpp HamiltonPath.cpp
UNIX> g++ Hamilton-Main.cpp 
UNIX> a.out 0
4
UNIX> a.out 1
0
UNIX> a.out 2
0
UNIX> a.out 3
24
UNIX>
And submit!