Add files via upload

NTT_fromR.nb

@@ -0,0 +1,89 @@
+(* Content-type: application/vnd.wolfram.mathematica *)
+
+(*** Wolfram Notebook File ***)
+(* http://www.wolfram.com/nb *)
+
+(* CreatedBy='Mathematica 11.2' *)
+
+(*CacheID: 234*)
+(* Internal cache information:
+NotebookFileLineBreakTest
+NotebookFileLineBreakTest
+NotebookDataPosition[       158,          7]
+NotebookDataLength[      2708,         81]
+NotebookOptionsPosition[      2399,         67]
+NotebookOutlinePosition[      2755,         83]
+CellTagsIndexPosition[      2712,         80]
+WindowFrame->Normal*)
+
+(* Beginning of Notebook Content *)
+Notebook[{
+Cell[BoxData[{
+ RowBox[{"SetDirectory", "[", 
+  "\"\</Users/maryam/Google Drive/Research/RF_improvment\>\"", 
+  "]"}], "\[IndentingNewLine]", 
+ RowBox[{"<<", "NumberTheoryFunctions.m"}], "\[IndentingNewLine]", 
+ RowBox[{"<<", " ", "numberth.m"}], "\[IndentingNewLine]", 
+ RowBox[{"Print", "[", "1", "]"}], "\[IndentingNewLine]", 
+ RowBox[{"init", "[", "]"}], "\n", 
+ RowBox[{
+  RowBox[{"U", " ", "=", " ", 
+   RowBox[{"ReadList", "[", "\n", "   ", 
+    RowBox[{
+    "\"\</Users/maryam/Google Drive/Research/RF_improvment/testNTT.txt\>\"", 
+     ",", " ", "\n", "   ", "Number", ",", " ", 
+     RowBox[{"RecordLists", " ", "->", " ", "True"}]}], "]"}]}], 
+  ";"}], "\n", 
+ RowBox[{
+  RowBox[{"A", " ", "=", " ", 
+   RowBox[{"U", "[", 
+    RowBox[{"[", "1", "]"}], "]"}]}], ";"}], "\n", 
+ RowBox[{
+  RowBox[{"B", " ", "=", " ", 
+   RowBox[{"U", "[", 
+    RowBox[{"[", "2", "]"}], "]"}]}], ";"}], "\n", 
+ RowBox[{
+  RowBox[{"t1", " ", "=", " ", 
+   RowBox[{
+    RowBox[{"fnt", "[", "A", "]"}], "*", 
+    RowBox[{"fnt", "[", "B", "]"}]}]}], ";"}], "\n", 
+ RowBox[{
+  RowBox[{"c1", " ", "=", " ", 
+   RowBox[{"ifnt", "[", "t1", "]"}]}], ";"}], "\n", "c1", "\n", 
+ RowBox[{"str", " ", "=", " ", 
+  RowBox[{
+  "OpenWrite", "[", "\n", "  ", 
+   "\"\</Users/maryam/Google Drive/Research/RF_improvment/outNTT.txt\>\"", 
+   "]"}]}], "\n", 
+ RowBox[{"Write", "[", 
+  RowBox[{"str", ",", " ", "c1"}], "]"}], "\n", 
+ RowBox[{"Close", "[", "str", "]"}]}], "Input",
+ CellChangeTimes->{{3.735253411577821*^9, 3.735253456196847*^9}, {
+   3.735253572224071*^9, 3.73525357352667*^9}, 3.7352543189796953`*^9, {
+   3.735267261221972*^9, 3.735267263104089*^9}, {3.735268133707547*^9, 
+   3.735268166321697*^9}, 3.735268365769558*^9, 3.735268476400447*^9, 
+   3.735268511237224*^9},
+ CellLabel->"In[32]:=",ExpressionUUID->"f446b4b9-745a-4d9f-8637-212fe444c6d8"]
+},
+WindowSize->{808, 713},
+WindowMargins->{{Automatic, 297}, {-46, Automatic}},
+FrontEndVersion->"11.3 for Mac OS X x86 (32-bit, 64-bit Kernel) (March 5, \
+2018)",
+StyleDefinitions->"Default.nb"
+]
+(* End of Notebook Content *)
+
+(* Internal cache information *)
+(*CellTagsOutline
+CellTagsIndex->{}
+*)
+(*CellTagsIndex
+CellTagsIndex->{}
+*)
+(*NotebookFileOutline
+Notebook[{
+Cell[558, 20, 1837, 45, 367, "Input",ExpressionUUID->"f446b4b9-745a-4d9f-8637-212fe444c6d8"]
+}
+]
+*)
+

NTT_fromR.wls

@@ -0,0 +1,17 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+SetDirectory["/Users/maryam/Google Drive/Research/RF_improvment"]
+<<NumberTheoryFunctions.m
+<< numberth.m
+
+init[]
+U = ReadList[
+   "/Users/maryam/Google Drive/Research/RF_improvment/testNTT.txt", 
+   Number, RecordLists -> True];
+A = U[[1]];
+B = U[[2]];
+t1 = fnt[A]*fnt[B];
+c1 = ifnt[t1];
+c1;
+Export["/Users/maryam/Google Drive/Research/RF_improvment/outNTT.txt",c1]

RF_FFT.R

@@ -0,0 +1,126 @@
+library("ape")
+library("phangorn")
+library("readtext")
+
+#Beta function
+Beta=function(m){
+  if(m<0){return(1)}
+  return(dfactorial(2*m+1))
+}
+
+# the function for computing the number of internal edges
+
+internaledges <- function(tree,ntip){
+  intedges=array(0,c(1,ntip-1))
+  edges=tree$edge
+  for (i in (2*ntip-1):(ntip+1)) {
+    children=which(edges[,1]==i)
+    child1=edges[children[1],2]
+    child2=edges[children[2],2]
+    if((child1 <= ntip)&(child2 <= ntip)){intedges[i-ntip]=0}
+    else if((child1<= ntip) & (child2 > ntip)) {intedges[i-ntip]=intedges[child2-ntip]+1}
+    else if((child2<= ntip )& (child1 > ntip)){intedges[i-ntip]=intedges[child1-ntip]+1}
+    else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2}
+  }
+  return(intedges)
+}
+
+
+internalchildren <- function(tree,v,ntip){
+  edges=tree$edge
+  children=which(edges[,1]==v)
+  child1=edges[children[1],2]
+  child2=edges[children[2],2]
+  if((child1 > ntip) & (child2 > ntip)){result=c(2,child1,child2)}
+  else if((child1 > ntip) & (child2 <= ntip)){result=c(1,child1)}
+  else if((child2 > ntip) & (child1 <= ntip)){result=c(1,child2)}
+  else {result=0}
+  return(result)
+}
+
+
+RF_Convolve=function(tree,n){
+  ntip=n-1
+  N=tree$Nnode
+  R=rep(list(matrix(0,(ntip-1),(ntip-1))),N)
+  edges=internaledges(tree,ntip)
+  B=c()
+  for (k in 0:(n-2)) {
+    B[k+1]=Beta(k)        
+  }
+  for (v in N:1) {
+    intchild=internalchildren(tree,v+ntip,ntip)
+    intedges=edges[v]
+    if(intchild[1]==0){
+      R[[v]][1,1]=1
+    }
+    else if(intchild[1]==1){
+      Rchild=R[[intchild[2]-ntip]]
+      R[[v]][1,intedges+1]=1
+      R[[v]][2:(ntip-1),1]=rowSums(t(t(Rchild[1:(ntip-2),])*B[1:(ntip-1)]))
+      R[[v]][2:(ntip-1),2:(ntip-1)]=Rchild[2:(ntip-1),1:((ntip-2))]
+    }
+    else {
+      Rchild1=R[[intchild[2]-ntip]]
+      Rchild2=R[[intchild[3]-ntip]]
+      R[[v]][1,intedges+1]=1
+      R[[v]][3,1]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))
+      for (s in 4:(ntip-1)) {
+        R[[v]][s,1]=sum(rowSums(t(t(Rchild1[1:(s-2),])*B[1:(ntip-1)]))*rowSums(t(t(Rchild2[(s-2):1,])*B[1:(ntip-1)])))
+      }
+      sum1=matrix(0,(ntip-2),(ntip-2))
+      sum1[1,1:(ntip-2)]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*Rchild2[1,1:(ntip-2)]
+      for (s in 3:(ntip-1)) {
+        temp=colSums(rowSums(t(t(Rchild1[1:(s-1),])*B[1:(ntip-1)]))*Rchild2[(s-1):1,1:(ntip-2)])
+        sum1[s-1,1:(ntip-2)]=temp
+      }
+      sum2=matrix(0,(ntip-2),(ntip-2))
+      sum2[1,1:(ntip-2)]=sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))*Rchild1[1,1:(ntip-2)]
+      for (s in 3:(ntip-1)) {
+        temp=colSums(rowSums(t(t(Rchild2[1:(s-1),])*B[1:(ntip-1)]))*Rchild1[(s-1):1,1:(ntip-2)])
+        sum2[s-1,1:(ntip-2)]=temp
+      }
+
+      R1=Rchild1[1:(ntip-1),1:(ntip-3)]
+      R2=Rchild2[1:(ntip-1),1:(ntip-3)]
+      R1aug=cbind(R1,matrix(0,nrow(R1),ncol(R1)))
+      R2aug=cbind(R2,matrix(0,nrow(R2),ncol(R2)))
+      U=round(convolve(t(R1aug), rev(t(R2aug)), type="open"),0)
+      sum3=t(matrix(c(U,0),2*ncol(R1),2*nrow(R1))[1:ncol(R1),1:nrow(R1)])
+      sum3=cbind(array(0, dim=c(nrow(R1)-1,1)),sum3[2:nrow(R1),])
+      R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3
+    }
+  }
+  return(R)
+}
+
+
+
+#==========================================
+RsT=function(R,n,s){
+  B=c()
+  for (k in 0:(n-2)) {
+    B[k+1]=Beta(k)        
+  }
+  rst =sum(t(t(R[[1]][s+1,1:(n-2-s)])*B[1:(n-2-s)]))
+  return(rst)
+}
+
+#Compute the value of q_m(T)
+qmT=function(R,n,m){
+  qmt=0
+  for (s in m:(n-3)) {
+    rst=RsT(R,n,s)
+    qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m)
+  }
+  return(qmt)
+}
+
+polynomial=function(tree,n){
+  R=RF_Convolve(tree,n)
+  for (i in seq(0,2*(n-3),2)) {
+    print(qmT(R,n,n-3-(i/2)))
+  }
+}
+
+

RF_NTT.R

@@ -0,0 +1,142 @@
+library("ape")
+library("phangorn")
+library("readtext")
+
+#Beta function
+Beta=function(m){
+  if(m<0){return(1)}
+  return(dfactorial(2*m+1))
+}
+
+# the function for computing the number of internal edges
+internaledges <- function(tree,ntip){
+  intedges=array(0,c(1,ntip-1))
+  edges=tree$edge
+  for (i in (2*ntip-1):(ntip+1)) {
+    children=which(edges[,1]==i)
+    child1=edges[children[1],2]
+    child2=edges[children[2],2]
+    if(child1 <= ntip&child2 <= ntip){intedges[i-ntip]=0}
+    else if(child1<= ntip & child2 > ntip){intedges[i-ntip]=intedges[child2-ntip]+1}
+    else if(child2<= ntip & child1 > ntip){intedges[i-ntip]=intedges[child1-ntip]+1}
+    else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2}
+  }
+  return(intedges)
+}
+
+# the function for computing the number of internal children
+internalchildren <- function(tree,v,ntip){
+  edges=tree$edge
+  children=which(edges[,1]==v)
+  child1=edges[children[1],2]
+  child2=edges[children[2],2]
+  if(child1 > ntip & child2 > ntip){result=c(2,child1,child2)}
+  else if(child1 > ntip & child2 <= ntip){result=c(1,child1)}
+  else if(child2 > ntip & child1 <= ntip){result=c(1,child2)}
+  else {result=0}
+  return(result)
+}
+
+RF_Convolve=function(tree,n){
+  ntip=n-1
+  N=tree$Nnode
+  R=rep(list(matrix(0,(ntip-1),(ntip-1))),N)
+  edges=internaledges(tree,ntip)
+  B=c()
+  for (k in 0:(n-2)) {
+    B[k+1]=Beta(k)        
+  }
+  for (v in N:1) {
+    intchild=internalchildren(tree,v+ntip,ntip)
+    intedges=edges[v]
+    if(intchild[1]==0){
+      R[[v]][1,1]=1
+    }
+    else if(intchild[1]==1){
+      Rchild=R[[intchild[2]-ntip]]
+      R[[v]][1,intedges+1]=1
+      R[[v]][2:(ntip-1),1]=rowSums(t(t(Rchild[1:(ntip-2),])*B[1:(ntip-1)]))
+      R[[v]][2:(ntip-1),2:(ntip-1)]=Rchild[2:(ntip-1),1:((ntip-2))]
+    }
+    else {
+      Rchild1=R[[intchild[2]-ntip]]
+      Rchild2=R[[intchild[3]-ntip]]
+      R[[v]][1,intedges+1]=1
+      R[[v]][3,1]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))
+      for (s in 4:(ntip-1)) {
+        R[[v]][s,1]=sum(rowSums(t(t(Rchild1[1:(s-2),])*B[1:(ntip-1)]))*rowSums(t(t(Rchild2[(s-2):1,])*B[1:(ntip-1)])))
+      }
+      sum1=matrix(0,(ntip-2),(ntip-2))
+      sum1[1,1:(ntip-2)]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*Rchild2[1,1:(ntip-2)]
+      for (s in 3:(ntip-1)) {
+        temp=colSums(rowSums(t(t(Rchild1[1:(s-1),])*B[1:(ntip-1)]))*Rchild2[(s-1):1,1:(ntip-2)])
+        sum1[s-1,1:(ntip-2)]=temp
+      }
+      sum2=matrix(0,(ntip-2),(ntip-2))
+      sum2[1,1:(ntip-2)]=sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))*Rchild1[1,1:(ntip-2)]
+      for (s in 3:(ntip-1)) {
+        temp=colSums(rowSums(t(t(Rchild2[1:(s-1),])*B[1:(ntip-1)]))*Rchild1[(s-1):1,1:(ntip-2)])
+        sum2[s-1,1:(ntip-2)]=temp
+      }
+
+      R1=Rchild1[1:(ntip-1),1:(ntip-3)]
+      #R1aug=numeric(nrow(R1)*(2*ncol(R1)-1))
+      R1aug=numeric(2048)
+      t=1
+      for(i in 1:ncol(R1)){
+        R1aug[t:(t+nrow(R1)-1)]=R1[,i]
+        t=t+3*nrow(R1)
+      }
+      R2=Rchild2[1:(ntip-1),1:(ntip-3)]
+      #R2aug=numeric(nrow(R2)*(2*ncol(R2)-1))
+      R2aug=numeric(2048)
+      t=1
+      for(i in 1:ncol(R2)){
+        R2aug[t:(t+nrow(R2)-1)]=R2[,i]
+        t=t+3*nrow(R2)
+      }
+
+      write(R1aug,"testNTT.txt",ncolumns=2048,append = TRUE)
+      write(R2aug,"testNTT.txt",ncolumns=2048,append = TRUE)
+     #read the output of mathematica
+      system("/Applications/Mathematica.app/Contents/MacOS/MathKernel -script  /Users/maryam/Desktop/NTT_fromR.wls")
+      U=as.matrix(read.csv("/Users/maryam/Google Drive/Research/RF_improvment/outNTT.txt",header = FALSE, sep="'", quote=""))
+      Matc=ceiling(length(R1aug)/(3*nrow(R1)))
+      sum3=matrix(c(U,numeric(Matc*3*nrow(R1)-length(R1aug))),nrow=3*nrow(R1))[1:nrow(R1),1:ncol(R1)]
+      sum3=cbind(array(0, dim=c(nrow(R1)-1,1)),sum3[2:nrow(R1),])
+      R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3
+      file.remove("testNTT.txt")
+    }
+  }
+
+  return(R)
+}
+
+
+
+#==========================================
+RsT=function(R,n,s){
+  B=c()
+  for (k in 0:(n-2)) {
+    B[k+1]=Beta(k)        
+  }
+  rst =sum(t(t(R[[1]][s+1,1:(n-2-s)])*B[1:(n-2-s)]))
+  return(rst)
+}
+
+#Compute the value of q_m(T)
+qmT=function(R,n,m){
+  qmt=0
+  for (s in m:(n-3)) {
+    rst=RsT(R,n,s)
+    qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m)
+  }
+  return(qmt)
+}
+
+polynomial=function(tree,n){
+  R=RF_Convolve(tree,n)
+  for (i in seq(0,2*(n-3),2)) {
+    print(qmT(R,n,n-3-(i/2)))
+  }
+}