Some progress on the Hilbert series of CY orbifolds.
(1) The general formula is

PI + sum_{B in BB} Porb(B,0) + sum_{C in CC} PC(C)

here (1) the **initial term** is

PI = (1 + at + bt^2 + at^3 + t^4)/(1-t)^4

in the point sum, each B = 1/r(a1,a2,a3), and the term is
  Porb(B,0) = Num/(1-t)^3/(1-t^r)
where, in the isolated case, Num is the unique polynomial with support in
[3,r] which is the inverse of prod ((1-t^ai)/(1-t)); in general you also
have to take out /(1-t^b) for common factors, etc., and the function is
given by Magma function Qorb below.

in the curve sum, each C = 1/r(a,r-a) (+ data), and the term is
  PC = Num/(1-t)^2(1-t^r)^2
where Num is a symmetric polynomial of degree 2r+2 supported in [3,2r-1].
Num has r-1 arbitrary coefficients, and my current working assumption is
that all values in an open range will occur.

(2) I have a routine to reverse engineer this formula for each Y in the
list of 7555 qsmooth hypersurface. Here P(A) and PI(A) are obvious things.
Bask(A) is selfcontained; it just calculates the list [ai mod r, d mod r]
for A = [a1,..a5,d], and puts the singular points and curves it finds
into baskets B and C. PointTerms is the sum over point terms, using the
preprogrammed function Qorb. X(A) is the sum of the curve terms, treated
here as an unknown to be taken apart. X(A) is supposed to be a sum of
terms Num_r/Denom([1,1,r,r]) over the baskets. The function PC determines
the different numerators Num_r using PartialFractionDecomposition plus
a couple of little twists -- this is recent and only tested in a few
hundred cases. It may go wrong, for example if Y has curves of
singularities with index having a common factor.

See below for examples of its use.

// Given a qsmooth CY hypersurface Y(d) in PP(a1,.. a5)
// in the form A := [a1,..a5,d], parse it into the form
// PI (initial term)
// plus sum Porb(1/r(a,b,c)) with 1/r(a,b,c) in PointBasket
// plus sum Ar/[1,1,r,r] where Ar is a symmetric polynomial
// of symmetric degree 2r+2 with support in [3..2r-1].

function P(A)
return (1-t^A[6])/Denom(A[1..5]);
end function;

function PI(A) // Initial term P_I, only for CY 3folds
n1 := #[i : i in A | i eq 1];
n2 := #[i : i in A | i eq 2];
return (1 + (n1-4)*t + (n2+Binomial(n1-3,2))*t^2 + (n1-4)*t^3 + t^4)/(1-t)^4;
end function;

function Bask(A);
 B := []; C := [];
 Relevant := [r : r in [2..A[5]] | r eq GCD([Integers()
             | A[i] : i in [1..5] | A[i] mod r eq 0])];
 for r in Relevant do
  Amod := [a mod r : a in A ];
  case [#[a : a in Amod[1..5] | a eq 0],#[a : a in Amod[6..6] | a eq 0]]:
   // when 0:; // no sing, do nothing when 4:; //error when 5:; //error
   when [1,0]:
    Append(~B, Insert(Exclude(Exclude(Amod[1..5],Amod[6]),0),1,r));
    // end this case: If ai divides d then Pi not on X
   when [2,1]:
     // calculate the number S^0(d) in PP(ai,aj)
     Num := Floor(A[6]/LCM([a : a in A[1..5] | a mod r eq 0]));
    for i in [1..Num] do
    Append(~B, Insert(Exclude(Exclude(Amod[1..5],0),0),1,r));
    end for;
    // end if; If r does not divides d then Lij is line of 1/r
   when [2,0]:
    Append(~C, Insert(Exclude(Exclude(Exclude(Amod[1..5],0),0),Amod[6]),1,r));
   when [3,1]:; // necessarily r divides d and curve of 1/r
    Append(~C, Insert(Exclude(Exclude(Exclude(Amod[1..5],0),0),0),1,r));
  end case;
end for;
return B, C;
end function;

function PointTerms(B)
return &+[K | Qorb(b[1],b[2..4],0) : b in B];
end function;

function X(A)
return P(A)-PI(A)-PointTerms(Bask(A));
end function;

function PC(A) // The curve terms P_C
B,C := Bask(A);
YY := PartialFractionDecomposition(X(A)/t^3*(1-t)^4);
return [t^3/(1-t)^4*&+[K|y[3]/y[1]^y[2] : y in YY
  | IsDivisibleBy(1-t^r,y[1])] where r is c[1] : c in C];
end function;

===================

for example, a little segment from the famous 7555:

AA:=[[1,12,27,32,36,108],[1,12,27,40,40,120],[1,12,27,68,96,204],[1,12,27,68,108,216],[1,12,32,39,45,129],[1,12,33,40,46,132],[1,12,33,92,138,276],[1,12,39,52,52,156],[1,12,39,52,65,169],[1,12,39,52,92,196],[1,12,39,52,103,207],[1,12,39,104,156,312],[1,12,40,93,134,280],[1,12,41,96,138,288],[1,12,42,98,141,294],[1,12,51,64,127,255],[1,12,51,64,128,256],[1,12,53,120,186,372],[1,12,54,68,81,216],[1,12,54,122,189,378],[1,12,64,153,230,460],[1,12,66,92,105,276],[1,12,66,158,237,474],[1,12,76,177,266,532],[1,12,77,180,270,540],[1,12,78,104,117,312],[1,13,15,30,46,105],[1,13,21,28,62,125],[1,13,21,35,35,105],[1,13,21,35,69,139],[1,13,22,29,52,117],[1,13,23,28,32,97],[1,13,28,56,98,196],[1,13,28,70,111,223],[1,13,29,43,73,159],[1,13,29,44,74,161],[1,13,30,75,106,225],[1,13,34,41,88,177],[1,13,34,61,108,217],[1,13,34,95,142,285],[1,13,41,96,150,301],[1,13,41,109,163,327]];

> AA[34];
[ 1, 13, 28, 70, 111, 223 ]
> Bask(AA[34]);
[
    [ 13, 1, 5, 7 ],
    [ 28, 1, 13, 14 ],
    [ 70, 1, 28, 41 ],
    [ 111, 13, 28, 70 ]
]
[
    [ 14, 1, 13 ]
]
> X(AA[34]);
(-t^22 - t^20 - t^18 - t^17 - t^16 - t^15 - t^14 - t^13 - t^12 - t^10 -
    t^8)/(t^30 - 2*t^29 + t^28 - 2*t^16 + 4*t^15 - 2*t^14 + t^2 - 2*t + 1)
> PartialFractionDecomposition(X(AA[34])/t^3*(1-t)^4);
[
    <t + 1, 2, 5/49>,
    <t^6 - t^5 + t^4 - t^3 + t^2 - t + 1, 1, 29/196*t^4 - 87/196*t^3 +
        125/196*t^2 - 47/98*t - 3/98>,
    <t^6 - t^5 + t^4 - t^3 + t^2 - t + 1, 2, 1/7*t^5 - 5/28*t^3 + 1/7*t^2 -
        1/7*t + 5/28>,
    <t^6 + t^5 + t^4 + t^3 + t^2 + t + 1, 1, -1/4*t^4 + 1/4*t^3 - 1/4*t^2>,
    <t^6 + t^5 + t^4 + t^3 + t^2 + t + 1, 2, -1/2*t^5 - 1/4*t^3 - 1/4>
]

> [X(A) eq &+[K| x : x in PC(A)] : A in AA];
// gives "true" 42 times, and checks that the routines work without crashing.

===================

(III) The above analyses existing examples, and can be quickly programmed
to put the 7555 hypersurfaces into a database. I want to use these ideas
to predict new example. We expect the nicely generated CYs to be sparse
in the set of all possible Hilbert series, so that if I modify P(Y) by a
clumsy amount, it is unlikely that the new P(Y) will correspond to a Y
that we can work with. However, I can modify the above Num_r to make more
delicate variations, e.g.,

> A := [1,4,5,5,5,20];
is a hypersurface with a curve Ga of 1/5(1,4), and Ga contributes
  (2*t^9 + 4*t^8 + 6*t^7 + 2*t^6 + 6*t^5 + 4*t^4 + 2*t^3)/[1,1,5,5]
I modify it by changing the 6,2,6 in the middle to 5,4,5 (which should
remove a generator in deg 5, leave the point singularities and the
degree of the curve unchanged).

> P1 := (1-t^20)/Denom([1,4,5,5,5]);
  // the Hilbert series of Y(20) in PP(1,4,5,5,5)
> P2 := P1 + (-t^7+2*t^6-t^5)/Denom([1,1,5,5]);
> P2*Denom([1,4,5,5,6,9]);
  t^30 - t^18 - t^12 + 1

That is, the modification gives Y(12,18) in PP(1,4,5,5,6,9)

Or change the 4,6,2,6,4 in the middle to 3,6,4,6,3
(which should remove a generator in deg 4 and add one in deg 6).
> P3 := P1 + (-t^8+2*t^6-t^4)/Denom([1,1,5,5]);
gives a plausible codim 4 guy with 9 x 16 resolution

Or change the Numerator of the Ga term to
  (2*t^9 + 3*t^8 + 7*t^7 + 2*t^6 + 7*t^5 + 3*t^4 + 2*t^3)

> P4 := P1 - (t^8+-t^7-t^5+t^4)/Denom([1,1,5,5]);
> P4*Denom([1,5,5,5,7,8,9]);
  -t^40 + t^26 + t^25 + t^24 + t^23 + t^22 - t^18 - t^17 - t^16 - t^15\
  - t^14 + 1
gives candidate codim 3 Y in PP(1,5,5,5,7,8,9) with 5x5 Pfaffian matrix of
degrees
    5,6,7,8
      7,8,9
        9,10
          11

We can play this game about a million times with the existing list.

===================
> load "CY/AAKS2";
// Loading "CY/AAKS2"
> #AAKS2;
// 14817

function BB(AA,r,a) // find transverse 1/r(a,r-a) curves in list AA of CYs
 BB0 := [A : A in AA | (([r,a,r-a] in C) or ([r,r-a,a] in C)) where B,C is Bask(A)];
 BB1 := [BB0[i] : i in [1..#BB0] | &and[GCD(b[1],r) eq 1 : b in Bask(BB0[i])]];
 return [BB1[i] : i in [1..#BB1] | (BB1[i][6] mod r) in [0,a,r-a]];
end function;

> AA := BB(AAKS2,5,2);
> #AA;
// 19

[ 4, 5, 6, 65, 66, 67, 69 ]