#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 2005;
int dp[MAX_N][2][MAX_N][3][2];
int sb[MAX_N];
vector<int> G[MAX_N];
int add(int x, int y) { return (x + y) % 1000000007; }
int sub(int x, int y) { return (x + 1000000007 - y) % 1000000007; }
int mul(int x, int y) { return (long long)x * y % 1000000007; }
int inv[MAX_N], fac[MAX_N], finv[MAX_N];
void make() {
  fac[0] = fac[1] = 1;
  finv[0] = finv[1] = 1;
  inv[1] = 1;
  for (int i = 2; i < MAX_N; i++) {
    inv[i] = 1000000007 -
             (long long)inv[1000000007 % i] * (1000000007 / i) % 1000000007;
    fac[i] = (long long)fac[i - 1] * i % 1000000007;
    finv[i] = (long long)finv[i - 1] * inv[i] % 1000000007;
  }
}
int comb(int a, int b) {
  if (a < b) {
    return 0;
  }
  return fac[a] * ((long long)finv[b] * finv[a - b] % 1000000007) % 1000000007;
}
int dfs(int u, int p) {
  sb[u]++;
  dp[u][0][1][2][0] = 1;
  int mx = 1;
  for (auto&(v) : (G[u])) {
    if (v != p) {
      sb[u] += dfs(v, u);
      for (int i = (int)(1); i < (int)(mx + 1); ++i) {
        for (int j = sb[v]; j >= 1; j--) {
          for (int l = 0; l < (int)(2); ++l) {
            for (int m = 0; m < (int)(2); ++m) {
              dp[u][1][i + j][0][l ^ m] =
                  add(dp[u][1][i + j][0][l ^ m],
                      mul(dp[u][0][i][0][l],
                          mul(add(dp[v][0][j][0][m], dp[v][0][j][1][m]) +
                                  dp[v][0][j][2][m],
                              comb(i + j, i))));
            }
          }
        }
        for (int j = sb[v]; j >= 1; j--) {
          for (int l = 0; l < (int)(2); ++l) {
            for (int m = 0; m < (int)(2); ++m) {
              dp[u][1][i + j - 1][0][l ^ m ^ 1] =
                  add(dp[u][1][i + j - 1][0][l ^ m ^ 1],
                      mul(dp[u][0][i][1][l],
                          mul(dp[v][0][j][2][m], comb(i + j - 1, i))));
              dp[u][1][i + j - 1][0][l ^ m ^ 1] = add(
                  dp[u][1][i + j - 1][0][l ^ m ^ 1],
                  mul(dp[u][0][i][1][l],
                      mul(mul(dp[v][0][j][1][m], inv[2]), comb(i + j - 1, i))));
            }
          }
        }
        for (int j = sb[v]; j >= 1; j--) {
          for (int l = 0; l < (int)(2); ++l) {
            for (int m = 0; m < (int)(2); ++m) {
              dp[u][1][i + j][1][l ^ m] =
                  add(dp[u][1][i + j][1][l ^ m],
                      mul(dp[u][0][i][1][l],
                          mul(dp[v][0][j][0][m], comb(i + j, i))));
              dp[u][1][i + j][1][l ^ m] =
                  add(dp[u][1][i + j][1][l ^ m],
                      mul(dp[u][0][i][1][l],
                          mul(mul(mul(dp[v][0][j][1][m], inv[2]), i),
                              comb(i + j, i + 1))));
              dp[u][1][i + j][1][l ^ m] =
                  add(dp[u][1][i + j][1][l ^ m],
                      mul(dp[u][0][i][1][l],
                          mul(mul(dp[v][0][j][1][m], inv[2]), comb(i + j, i))));
              dp[u][1][i + j][1][l ^ m] =
                  add(dp[u][1][i + j][1][l ^ m],
                      mul(dp[u][0][i][1][l],
                          mul(mul(dp[v][0][j][2][m], i), comb(i + j, i + 1))));
            }
          }
        }
        for (int j = sb[v]; j >= 1; j--) {
          for (int l = 0; l < (int)(2); ++l) {
            for (int m = 0; m < (int)(2); ++m) {
              dp[u][1][i + j - 1][1][l ^ m ^ 1] =
                  add(dp[u][1][i + j - 1][1][l ^ m ^ 1],
                      mul(dp[u][0][i][2][l],
                          mul(2 * dp[v][0][j][2][m], comb(i + j - 1, i))));
              dp[u][1][i + j - 1][1][l ^ m ^ 1] =
                  add(dp[u][1][i + j - 1][1][l ^ m ^ 1],
                      mul(dp[u][0][i][2][l],
                          mul(dp[v][0][j][1][m], comb(i + j - 1, i))));
            }
          }
        }
        for (int j = sb[v]; j >= 1; j--) {
          for (int l = 0; l < (int)(2); ++l) {
            for (int m = 0; m < (int)(2); ++m) {
              dp[u][1][i + j][2][l ^ m] =
                  add(dp[u][1][i + j][2][l ^ m],
                      mul(dp[u][0][i][2][l],
                          mul(dp[v][0][j][0][m], comb(i + j, i))));
              dp[u][1][i + j][2][l ^ m] =
                  add(dp[u][1][i + j][2][l ^ m],
                      mul(dp[u][0][i][2][l],
                          mul(mul(dp[v][0][j][1][m], i), comb(i + j, i + 1))));
              dp[u][1][i + j][2][l ^ m] = add(
                  dp[u][1][i + j][2][l ^ m],
                  mul(dp[u][0][i][2][l],
                      mul(mul(dp[v][0][j][2][m], i - 1), comb(i + j, i + 1))));
            }
          }
        }
      }
      swap(dp[u][0], dp[u][1]);
      for (int j = 0; j < (int)(mx + 1); ++j) {
        for (int k = 0; k < (int)(2); ++k) {
          for (int l = 0; l < (int)(2); ++l) {
            dp[u][1][j][k][l] = 0;
          }
        }
      }
      mx += sb[v];
    }
  }
  return sb[u];
}
int main() {
  cin.tie(0);
  ios::sync_with_stdio(false);
  int n;
  cin >> n;
  for (int i = 0; i < (int)(n - 1); ++i) {
    int a;
    cin >> a;
    G[a - 1].push_back(i + 1);
  }
  make();
  dfs(0, -1);
  int ans = 0;
  for (int i = 0; i < (int)(n); ++i) {
    for (int j = 0; j < (int)(3); ++j) {
      ans = add(ans, dp[0][0][i + 1][j][0]);
      ans = sub(ans, dp[0][0][i + 1][j][1]);
    }
  }
  cout << ans << "\n";
  return 0;
}
