===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/tk_fd.rr,v
retrieving revision 1.7
retrieving revision 1.11
diff -u -p -r1.7 -r1.11
--- OpenXM/src/asir-contrib/packages/src/tk_fd.rr	2014/08/20 06:53:22	1.7
+++ OpenXM/src/asir-contrib/packages/src/tk_fd.rr	2014/08/27 02:14:40	1.11
@@ -1,8 +1,8 @@
-/* $OpenXM: OpenXM/src/asir-contrib/packages/src/tk_fd.rr,v 1.6 2014/08/09 09:27:48 takayama Exp $ */
+/* $OpenXM: OpenXM/src/asir-contrib/packages/src/tk_fd.rr,v 1.10 2014/08/26 06:59:56 takayama Exp $ */
 /*
   2015.05.23   h-mle/FD/Prog/fdpf.rr --> tk_fd.rr
  */
-#define USE_MODULE 0
+#define USE_MODULE 1
 import("names.rr")$
 import("ot_hgm_ahg.rr")$   // for check7
 import("ok_diff.rr")$  // for check7
@@ -188,6 +188,8 @@ localf hgpoly_fd1$
 localf fvect_poly2$
 localf check28$
 localf check28b$
+localf ahvec_beta$
+localf ahvec_abc$
 #endif
 /* Matsumoto, page 3 */
 
@@ -1969,6 +1971,7 @@ def ygdmat_abc(A,B,C) {
 }
 
 def ygfvect_poly(A,B,C,X) { 
+  if (C <=0 ) return(fvect_poly2(A,B,C,X));
   YgEvalDmat=1;
   if (type(X) == 4) X=newvect(length(X),X);
   if (taka_base_equal(B,YgFD_BB) && 
@@ -2014,7 +2017,7 @@ def ygahvec(A,B,C,Y) {
   Ye = 1;
   for (J=0; J<=M; J++) Ye = Ye*(Y[0][J]^E[0][J])*(Y[1][J]^E[1][J]);
   if (number_is_integer(A) && A < 0 && UseRecc) F=ygfvect_poly(A,B,C,X);
-  else F=ygfvect(A,B,C,X | option_list=Opt);  /* not implemented.*/
+  else F=fvect(A,B,C,X | option_list=Opt);  
   Fd=newvect(M+1);
   Fd[0] = F[0];  /* Fd are Euer derivatives of F_D */
   for (I=1; I<=M; I++) {
@@ -2035,12 +2038,12 @@ def ygahvec(A,B,C,Y) {
    /* Gamma factors, todo double check */
   if (NoGamma) {
      Gamma = "[1/gamma(e+1)]";
-  }else{
+  }else if (C>0) {
     Gamma=tk_number_invgamma(-A+1)*tk_number_invgamma(C-1+1);
     for (I=0; I<M; I++) {
       Gamma=Gamma*tk_number_invgamma(-B[I]+1);
     }
-  }
+  }else { Gamma=1; } /* ygfvect_poly2 is called. */
   R= [vtol(Ahvec),Gamma];
   if (All) return append(R,[Zvalue]);
   else return R;
@@ -2615,11 +2618,15 @@ def hgpoly_fd1(A,B,C) {
 def fvect_poly2(A,B,C,X) {
   if (C > 0) error("C is postive.");
   M=length(B);
-  CB = C-base_sum(B,0,0,length(B)-1)-1;
-  if (CB <= 0) error("c-sum(b)-1 is not positive.");
-  Steps = CB-1;
   if (type(X) == 4) X=newvect(length(X),X);
   Xrule = assoc(xvars(M),vtol(X));
+  CB = C-base_sum(B,0,0,length(B)-1)-1;
+  if (CB == 0) {
+    F = poly_fdvec(A,B,C)[0];
+    return(base_replace(F,Xrule));
+  }
+  if (CB < 0) error("c-sum(b)-1 is negative.");
+  Steps = CB-1;
   printf("Computing Dmat(ca) for parameters B=%a,X=%a\n",B,X);
   Dmat = ygDca2(a,B,c,X);
   Umat = map(red,matrix_inverse(Dmat));
@@ -2656,6 +2663,77 @@ def check28b(T) {
   return(F1-F2);
 }
 
+/*
+ works in case of row sum is zero or column sum is 0.
+*/
+def ahvec_beta(Beta,Y) {
+  Opt = getopt();
+  Ma=beta2marginal(Beta);
+  M=length(Beta)-2;
+  RS=newvect(2,Ma[0]);
+  CS=newvect(M+1,Ma[1]);
+  Mnew=M;
+  Fvec=newvect(M+1);
+  for (I=0, J=0; I<=M; I++) {
+    if (CS[I]==0) {
+      Mnew--; Fvec[I] = -1;
+    }else{ Fvec[I] = J; J++; }
+  }
+  if (Mnew < 0) error("zero matrix");
+  Ynew=newmat(2,Mnew+1);    
+  CSnew=newvect(Mnew+1);
+  for (I=0; I<=M; I++) {
+    if (Fvec[I] >= 0) {
+      J = Fvec[I];
+      Ynew[0][J] = Y[0][I];
+      Ynew[1][J] = Y[1][I];
+      CSnew[J] = CS[I];
+    }
+  }
+  if (RS[1] == 0) {
+    T=RS[0]; RS[0] = RS[1]; RS[1] = T;
+    for (J=0; J<=Mnew; J++) {
+      T=Ynew[0][J];
+      Ynew[0][J] = Ynew[1][J];
+      Ynew[1][J] = T;
+    }
+  }
+  ABCnew=marginal2abc(vtol(RS),vtol(CSnew));
+  printf("RS=%a, CSnew=%a, Ynew=%a\n",RS,CSnew,Ynew);
+  Tnew = ygahvec(ABCnew[0],ABCnew[1],ABCnew[2],Ynew | option_list=Opt);
+  Fnew =Tnew[0];
+  for (I=0; I<=M; I++) {
+    if (Fvec[I] == -1) Fvec[I] = 0;
+    else {
+      J=Fvec[I];
+      Fvec[I] = Fnew[J];
+    }
+  }
+  return cons(Fvec,cdr(Tnew));
+}
+/*
+ahvec_beta([3,1,2,0,3],[[1,1/2,1/3,1/5],[1,1,1,1]]|all=1);
+
+[2402] abc2marginal(-5,[-4,-3],-2);
+[[5,4],[2,4,3]]
+[2403] marginal2abc([4,5],[2,4,3]);  
+[-4,[-4,-3],-1]   c is again negative.
+
+2014.08.26
+ahvec_beta([10,1,2,3,4],[[3,1,1,1],[1,1/2,1/3,1/5]]|all=1);
+returns "devision by 0".  Reason?
+
+ahvec_abc(-6,[-2,-3],-4,[[1,1/2,1/3],[1,1,1]] | all=1);  
+ -> 2014.08.27 it works.
+Todo: call ahvec_abc in cgi/r-fd.rr
+*/
+
+def ahvec_abc(A,B,C,X) {
+  Opt=getopt();
+  Ma=abc2marginal(A,B,C);
+  Beta=cons(Ma[0][0],Ma[1]);
+  return ahvec_beta(Beta,X | option_list=Opt);  
+}
 #if  !USE_MODULE
 setparam(2)$
 // setparam_by_number(2,2)$