Testing Update

MatthewStuber · Mar 25, 2018 · 0864319 · 0864319
1 parent c0501c7
commit 0864319
Show file tree

Hide file tree

Showing 5 changed files with 139 additions and 97 deletions.
diff --git a/src/lib/Parametric_Contractor.jl b/src/lib/Parametric_Contractor.jl
@@ -29,7 +29,6 @@ function PI_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
                      hj::Function,h::Function,
                      opt::Vector{Any},Eflag::Bool,
                      Iflag::Bool,eDflag::Bool) where {T<:AbstractFloat}
-  println(" ---------- begin out place newton ------------- ")
   # unpacks option file
   kmax::Int64 = opt[1]
   etol::Float64 = opt[2]
@@ -63,9 +62,6 @@ function PI_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
     B = Y*H
     M = Y*J
   end
-println("iteration ($k)")
-println("preconditioned H: $B")
-println("preconditioned J: $M")
 
   for i=1:nx
     S1 = Interval(0.0)
@@ -77,25 +73,19 @@ println("preconditioned J: $M")
         S2 += M[i,j]*(X[j]-x_mid[j])
       end
     end
-    #println("M[$i,$i].lo: $(M[i,i].lo)")
-    #println("M[$i,$i].hi: $(M[i,i].hi)")
-    #println("M check: $(M[i,i].lo*M[i,i].hi > 0.0)")
     if M[i,i].lo*M[i,i].hi > 0.0
       N[i] = x_mid[i] - (B[i]+S1+S2)/M[i,i]
     else
-      Ntemp = N
+      Ntemp = copy(N)
       eD,N[i],Ntemp[i] = extProcess(N[i],X[i],M[i,i],S1,S2,B[i],rtol)
       if eD == 1
         eDflag = true
-        Xtemp = X
+        Xtemp = copy(X)
         Xtemp[i] = Ntemp[i] ∩ X[i]
         X[i] = N[i] ∩ X[i]
         return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh
       end
     end
-    #println("N[$i]: $(N[i])")
-    #println("X[$i]: $(X[i])")
-    #println("N[$i] in X[$i]: $(Strict_XinY(N[i],X[i]))")
     if Strict_XinY(N[i],X[i])
       inclusion[i] = true
       inclusionHigh[i] = true
@@ -110,7 +100,6 @@ println("preconditioned J: $M")
         inclusionHigh[i] = true
       end
     end
-    #println("N[$i] & X[$i] not disjoint: $(~isdisjoint(N[i],X[i]))")
     if ~isdisjoint(N[i],X[i])
       X[i] = N[i] ∩ X[i]
     else
@@ -148,9 +137,6 @@ println("preconditioned J: $M")
       B = Y*H
       M = Y*J
     end
-  println("iteration ($k)")
-  println("preconditioned H: $B")
-  println("preconditioned J: $M")
 
     for i=1:nx
       S1 = Interval(0.0)
@@ -162,25 +148,19 @@ println("preconditioned J: $M")
           S2 += M[i,j]*(X[j]-x_mid[j])
         end
       end
-      #println("M[$i,$i].lo: $(M[i,i].lo)")
-      #println("M[$i,$i].hi: $(M[i,i].hi)")
-      #println("M check: $(M[i,i].lo*M[i,i].hi > 0.0)")
       if M[i,i].lo*M[i,i].hi > 0.0
         N[i] = x_mid[i] - (B[i]+S1+S2)/M[i,i]
       else
-        Ntemp = N
+        Ntemp = copy(N)
         eD,N[i],Ntemp[i] = extProcess(N[i],X[i],M[i,i],S1,S2,B[i],rtol)
         if eD == 1
           eDflag = true
-          Xtemp = X
+          Xtemp = copy(X)
           Xtemp[i] = Ntemp[i] ∩ X[i]
           X[i] = N[i] ∩ X[i]
           return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh
         end
       end
-      #println("N[$i]: $(N[i])")
-      #println("X[$i]: $(X[i])")
-      #println("N[$i] in X[$i]: $(Strict_XinY(N[i],X[i]))")
       if Strict_XinY(N[i],X[i])
         inclusion[i] = true
         inclusionHigh[i] = true
@@ -195,7 +175,6 @@ println("preconditioned J: $M")
           inclusionHigh[i] = true
         end
       end
-      #println("N[$i] & X[$i] not disjoint: $(~isdisjoint(N[i],X[i]))")
       if ~isdisjoint(N[i],X[i])
         X[i] = N[i] ∩ X[i]
       else
@@ -223,7 +202,6 @@ println("preconditioned J: $M")
     Eflag = true
   end
   Xtemp = copy(X)
-  println(" ---------- end out place newton ------------- ")
   return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh
 end
 
@@ -288,9 +266,6 @@ function PI_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
     B = Y*H
     M = eye(nx)-Y*J
   end
-  #println("iteration ($k)")
-  #println("preconditioned H: $B")
-  #println("preconditioned J: $M")
 
   for i=1:nx
     S1 = Interval(0.0)
@@ -353,9 +328,6 @@ function PI_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
       B = Y*H
       M = eye(nx)-Y*J
     end
-#    println("iteration ($k)")
-#    println("preconditioned H: $B")
-#    println("preconditioned J: $M")
 
     for i=1:nx
       S1 = Interval(0.0)
@@ -408,7 +380,6 @@ function PI_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
   if exclusion
     Eflag = true
   end
-#  println(" ----------end out place kraw ------------- ")
   return X,Eflag,Iflag,inclusionLow,inclusionHigh
 end
 
@@ -488,9 +459,6 @@ function PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
   h!(H,x_mid,P)
   hj!(J,X,P)
   Sparse_Precondition!(H,J,mid.(J),SSto)
-  println("iteration ($k)")
-  println("preconditioned H: $H")
-  println("preconditioned J: $J")
   Mt = transpose(J)
   for i=1:nx
     S1 = Interval(0.0)
@@ -508,11 +476,11 @@ function PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
     if S3.lo*S3.hi > 0.0
       N[i] = x_mid[i] - (H[i]+S1+S2)/S3
     else
-      Ntemp = N
+      Ntemp = copy(N)
       eD,N[i],Ntemp[i] = extProcess(N[i],X[i],S3,S1,S2,H[i],rtol)
       if eD == 1
         eDflag = true
-        Xtemp = X
+        Xtemp = copy(X)
         Xtemp[i] = Ntemp[i] ∩ X[i]
         X[i] = N[i] ∩ X[i]
         return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh
@@ -562,9 +530,6 @@ function PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
     h!(H,x_mid,P)
     hj!(J,X,P)
     Sparse_Precondition!(H,J,mid.(J),SSto)
-    println("iteration ($k)")
-    println("preconditioned H: $H")
-    println("preconditioned J: $J")
     Mt = transpose(J)
 
     for i=1:nx
@@ -583,11 +548,11 @@ function PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},
       if S3.lo*S3.hi > 0.0
         N[i] = x_mid[i] - (H[i]+S1+S2)/S3
       else
-        Ntemp = N
+        Ntemp = copy(N)
         eD,N[i],Ntemp[i] = extProcess(N[i],X[i],S3,S1,S2,H[i],rtol)
         if eD == 1
           eDflag = true
-          Xtemp = X
+          Xtemp = copy(X)
           Xtemp[i] = Ntemp[i] ∩ X[i]
           X[i] = N[i] ∩ X[i]
           return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh
@@ -669,7 +634,6 @@ function PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
                      hj!::Function,h!::Function,
                      opt::Vector{Any},Eflag::Bool,Iflag::Bool) where {T<:AbstractFloat}
 
-#  println(" ---------- begin in place kraw ------------- ")
   # unpacks option file
   kmax::Int64 = opt[1]
   etol::Float64 = opt[2]
@@ -707,9 +671,6 @@ function PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
   h!(H,x_mid,P)
   hj!(J,X,P)
   Sparse_Precondition!(H,J,mid.(J),SSto)
-#  println("iteration ($k)")
-#  println("preconditioned H: $H")
-#  println("preconditioned J: $J")
   Mt = transpose(J)
   for i=1:nx
     for q=(Mt.colptr[i]):(Mt.colptr[i+1]-1)
@@ -762,9 +723,6 @@ function PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
     h!(H,x_mid,P)
     hj!(J,X,P)
     Sparse_Precondition!(H,J,mid.(J),SSto)
-#    println("iteration ($k)")
-#    println("preconditioned H: $H")
-#    println("preconditioned J: $J")
     Mt = transpose(J)
     N = [Interval(0.0) for i=1:nx]
     for i=1:nx
@@ -817,6 +775,5 @@ function PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},
   if exclusion
     Eflag = true
   end
-#  println(" ---------- begin out place kraw ------------- ")
   return X,Eflag,Iflag,inclusionLow,inclusionHigh
 end
diff --git a/src/lib/Parametric_Test.jl b/src/lib/Parametric_Test.jl
@@ -1,3 +1,4 @@
+#=
 """
     Miranda(h::Function,X::Vector{Interval{T}},P::Vector{Interval{T}})
 
@@ -18,8 +19,8 @@ function Miranda(h::Function,X::Vector{Interval{T}},P::Vector{Interval{T}}) wher
   i = -1
   return i
 end
-
-
+=#
+#=
 """
     MirandaExc
 
@@ -42,7 +43,8 @@ function MirandaExc(h,X,P,Eflag_in)
    end
    return Eflag
 end
-
+=#
+#=
 """
 --------------------------------------------------------------------------------
 Function: partialIncTop
@@ -84,7 +86,8 @@ function partialIncTop(h,X,P,PIflag,incHigh)
   end
   return k,PIflag
 end
-
+=#
+#=
 """
 --------------------------------------------------------------------------------
 Function: partialIncBot
@@ -126,7 +129,8 @@ function partialIncBot(h,X,P,PIflag,incLow)
   end
   return k,PIflag
 end
-
+=#
+#=
 """
 --------------------------------------------------------------------------------
 Function: spectralR
@@ -158,7 +162,8 @@ function spectralR(hj,X,P,rho)
   end
   return k,eigval
 end
-
+=#
+#=
 """
 --------------------------------------------------------------------------------
 Function: BoundaryTest
@@ -214,7 +219,7 @@ function BoundaryTest(h,hj,X0,X,P,opt,PIcert,PIflag,Iflag,Eflag,inclusionLow,inc
   end
   return PIcert,PIflag,Iflag,Eflag
 end
-
+=#
 #=
 function NewBoundaryTest(Xtemp,X,P,Iflag)
   nx = length(X)

diff --git a/src/lib/Parametric_Utility.jl b/src/lib/Parametric_Utility.jl
@@ -144,10 +144,10 @@ end
 
 Subfunction to generate output for extended division.
 """
-function extDivide(A::Interval{T},B::Interval{T},C::Interval{T}) where {T<:AbstractFloat}
+function extDivide(A::Interval{T}) where {T<:AbstractFloat}
   if ((A.lo == -0.0) && (A.hi == 0.0))
-    B = Interval(-Inf,Inf)
-    C = B
+    B::Interval{T} = Interval(-Inf,Inf)
+    C::Interval{T} = B
     return 0,B,C
   end
   if (A.lo == 0.0)
@@ -163,7 +163,6 @@ function extDivide(A::Interval{T},B::Interval{T},C::Interval{T}) where {T<:Abstr
     C = Interval(1.0/A.hi,Inf)
     return 3,B,C
   end
-  return -1,B,C
 end
 
 """
@@ -177,40 +176,45 @@ function extProcess(N::Interval{T},X::Interval{T},Mii::Interval{T},
                     S1::Interval{T},S2::Interval{T},B::Interval{T},rtol::Float64) where {T<:AbstractFloat}
   v = 1
   Ntemp::Interval{T} = copy(N)
-  IML::Interval{T} = Interval(0.0)
-  IMR::Interval{T} = Interval(0.0)
   M::Interval{T} = (B+S1+S2)+Interval(-rtol,rtol)
-  if (M.lo<=0 || M.hi>=0)
-    N = Interval(-Inf,Inf)
-    return 0, N, Ntemp
+  if (M.lo<=0 && M.hi>=0)
+    println("branch1")
+    return 0, Interval(-Inf,Inf), Ntemp
   end
   if (v == 1)
-    k,IML,IMR = extDivide(Mii,IML,IMR)
+    println("Mii: $Mii")
+    k,IML::Interval{T},IMR::Interval{T} = extDivide(Mii)
+    println("k: $k")
     if (k == 1)
-      N = mid(X)-M*IML
-      return 0, N, Ntemp
+      println("branch2")
+      return 0, (mid(X)-M*IML), Ntemp
     elseif (k == 2)
-      N = mid(X)-M*IMR
-      return 0, N, Ntemp
+      println("branch3")
+      return 0, (mid(X)-M*IMR), Ntemp
     elseif (k == 3)
       NR = mid(X)-M*IMR
       NL = mid(X)-M*IML
+      print("NR: $NR")
+      print("NL: $NL")
       if (~isdisjoint(NL,X) && isdisjoint(NR,X))
-        N = NL
-        return 0, N, Ntemp
+        println("branch4")
+        return 0, NL, Ntemp
       elseif (~isdisjoint(NR,X) && isdisjoint(NL,X))
-        N = NR
-        return 0, N, Ntemp
+        println("branch5")
+        return 0, NR, Ntemp
       elseif (~isdisjoint(NL,X) && ~isdisjoint(NR,X))
+        println("branch6")
         N = NL
         Ntemp = NR
-        return 1, N, Ntemp
+        print("NR: $NR")
+        print("NL: $NL")
+        return 1, NL, NR
       else
+        println("branch7")
         return -1, N, Ntemp
       end
     end
-  else
-    N = Interval(-Inf,Inf)
-    return 0, N, Ntemp
   end
+  println("branch8")
+  return 0, N, Ntemp
 end