kboxes.py

# -*- coding: utf-8 -*-
# kboxes.py (Kaucher p-boxes)
#
# Kaucher p-boxes 
#
# P-boxes (probability boxes) defined in terms of Kaucher generalized intervals
#
# @author: Scott Ferson 
# Copyright © 2025 Scott Ferson. All rights reserved

# We're not sure about the best formulation for normal and triangular p-boxes.
# The calculus for generalized p-boxes turn on the dependence assumptions, and 
# in sometimes subtle ways.  Check the basic convolutions before freaking out 
# about backcalculating the As. Maraschin et al. are handling rounding wrongly.
# Enrique suggests https://vscodium.com and https://www.cursor.com

import math
import scipy.stats as sps
import matplotlib.pyplot as plt
from kaucher import *

def ln(a):  # throw this hack away
    if a<=0 : return(1)
    else: return(math.log(a))
    
def left(x) : 
    if isinstance(x,Kaucher): return x[0]
    if isinstance(x,(int,float)): return x
    if isinstance(x,(list,tuple) ): return min(x)

def right(x) : 
    if isinstance(x,Kaucher): return x[1]
    if isinstance(x,(int,float)): return x
    if isinstance(x,(list,tuple) ): return max(x)

def continued(x, terms=5):
  rest = x
  p = [0,1]; q = [1,0]   
  for i in range(terms):
    whole, frac = divmod(rest, 1)
    an = int(whole)
    pn = an*p[-1]+p[-2]
    qn = an*q[-1]+q[-2]
    if math.isclose(x, pn/qn) or 10 <= pn + qn : break
    p += [pn]; q += [qn]
    rest = 1/frac
  return(pn,qn) 

# STEPS (below), the number of discretization levels, is hard coded.  That is, 
# although you can change the value of STEPS on the fly, doing so doesn't alter 
# the default number of levels used by the constructors such as I(), U(), E(), 
# which remains at whatever value it was when the constructors were defined.
# But you can alter the number by specifying it when calling a constructor.
# For example, a = U(0,1,21) yields a len(a) = 21, 
#        [k(0.0, 0.0),
#         k(0.05, 0.05),
#         k(0.1, 0.1),
#         ...
#         k(0.95, 0.95),
#         k(1.0, 1.0)]

STEPS = 100 # number of discretization levels in a p-box
bOt = 0.001
tOp = 0.999

II = [i/STEPS for i in range(STEPS)]
JJ = [i/STEPS for i in range(1,STEPS+1)]
III = [bOt] + [i/STEPS for i in range(1,STEPS)]
JJJ = [i/STEPS for i in range(1,STEPS)] + [tOp]

def U(m:Kaucher, M:Kaucher, n=STEPS) : # note that STEPS is NOT evaluated, but the original value is used
    m = Kaucher(m)
    M = Kaucher(M)
    ii = [i/n for i in range(n)]
    jj = [i/n for i in range(1,n+1)]
    return([k(m[0]+i*(M[0]-m[0]),m[1]+j*(M[1]-m[1])) for i,j in zip(ii,jj)])

def E(m, n=STEPS):  # use E(1/lambda) for the other parameterization
    m = Kaucher(m)
    III = [bOt] + [_/n for _ in range(1,n)]
    JJJ = [_/n for _ in range(1,n)] + [tOp]
    L = sps.expon.ppf(III,scale=left(m))
    R = sps.expon.ppf(JJJ,scale=right(m))
    return [k(L[_],R[_]) for _ in range(n)]

def N2(m,s,n=STEPS) : # it's not at all clear this is the right formulation
    III = [bOt] + [_/n for _ in range(1,n)]
    JJJ = [_/n for _ in range(1,n)] + [tOp]
    L1 = sps.norm.ppf(III,loc=left(m),scale=left(s))
    L2 = sps.norm.ppf(III,loc=left(m),scale=right(s))
    R1 = sps.norm.ppf(JJJ,loc=right(m),scale=left(s))
    R2 = sps.norm.ppf(JJJ,loc=right(m),scale=right(s))
    return [k(min(L1[_],L2[_]), max(R1[_],R2[_])) for _ in range(n)]

def N(m,s,n=STEPS) : # it's not at all clear this is the right formulation
    III = [bOt] + [_/n for _ in range(1,n)]
    JJJ = [_/n for _ in range(1,n)] + [tOp]
    L1 = sps.norm.ppf(III,loc=left(m),scale=left(s))
    L2 = sps.norm.ppf(III,loc=left(m),scale=right(s))
    R1 = sps.norm.ppf(JJJ,loc=right(m),scale=left(s))
    R2 = sps.norm.ppf(JJJ,loc=right(m),scale=right(s))
    return [k(min(L1[_],L2[_]), max(R1[_],R2[_])) for _ in range(n)]

def T2(m:Kaucher,mid:Kaucher,M:Kaucher, n=STEPS):
    def pr(p,pm,m,mid,M):
        if p<=pm : return m + (p * (M - m) * (mid - m))**0.5
        else :     return M - ((1.0 - p) * (M - m) * (M - mid))**0.5
    m = Kaucher(m)
    mid = Kaucher(mid)
    M = Kaucher(M)
    pp = [q/(n-1) for q in range(n)]
    pm = (mid-m) / (M-m)
    c = []
    for p in pp:
        c = c + [Kaucher(pr(p,pm[0],m[0],mid[0],M[0]),pr(p,pm[1],m[1],mid[1],M[1]))]        
    return(c)   # dividing pm[0] and pm[1] prolly doesn't make sense  

def T(m:Kaucher,mid:Kaucher,M:Kaucher, n=STEPS):
    def triangular(m : float, mid : float, M : float) :
        def pr(p,pm,m,mid,M):
            if p<=pm : return m + (p * (M - m) * (mid - m))**0.5
            else :     return M - ((1.0 - p) * (M - m) * (M - mid))**0.5
        pp = [q/(n-1) for q in range(n)]
        pm = (mid-m) / (M-m)
        c = []
        for p in pp:
            c = c + [Kaucher(pr(p,pm,m,mid,M),pr(p,pm,m,mid,M))]        
        return(c)
    m = Kaucher(m)
    mid = Kaucher(mid)
    M = Kaucher(M)
    # brute force strategy , ... but it's clearly WRONG for improper parameters
    tt = [Kaucher(math.inf,-math.inf)] * STEPS 
    for i in [0,1] :
        for j in [0,1] :
            for k in [0,1] :
                tt = Env(tt,triangular(m[i],mid[j],M[k]))
    return tt

def I(m:Kaucher, M=None, n=STEPS) : # this is basically the minmax function
    m = Kaucher(m)
    if M is None : M = m
    else : M = Kaucher(M)
    pp = [q/(n-1) for q in range(n)]
    return([k(m[0],M[1]) for p in pp])    # not sure this is what it should be!
    # Plot(I(k(1,3)))
    # Plot(I(k(10,3)))
    # Plot(I(k(1,4), k(10,3))) # Should this yield [1,3]?  It wouldn't seem so.
    # # But what should it be?
    # a = I(k(1,4), k(10,3))
    # b = I(k(1,4).intersection(k(10,3))); Plot(a,lwd=6); Replot(b)
    # b = I(k(1,4).env(k(10,3))); Plot(a,lwd=6); Replot(b)
    
Empty = [Kaucher(math.inf,-math.inf)] * STEPS # degenerate vacuous distribution?

def Mid(a) :
    return (a[0].mid() + a[-1].mid()) / 2

def Straddles(a) :                                           # is this correct?
    L = a[0].range();  R = a[-1].range()
    LR = [L.a,L.b,R.a,R.b]
    return min(LR) <= 0 <= max(LR)
    
def Plot(a,fmt='',lab='',lwd='', PLT=None):
    if PLT is None : PLT = plt
    n = len(a)
    if fmt == '' : fmt = ('r-','k-')
    if isinstance(fmt, str) : fmt = (fmt,fmt)
    if lwd == '' : lwd = 1
    if not lab == '' : PLT.text(Mid(a),1.008,lab)
    for i in range(n-1) :
        x = [a[i][0],a[i][0],a[i+1][0]]
        p = [i/n,(i+1)/n,(i+1)/n]
        PLT.plot(x,p,fmt[0],lw=lwd) 
    for i in range(n-1) :
        x = [a[i][1],a[i+1][1],a[i+1][1]]
        p = [(i+1)/n,(i+1)/n,(i+2)/n]
        PLT.plot(x,p,fmt[1],lw=lwd)
    f = a[-1].is_proper()
    if f : PLT.plot([a[-1][0],a[-1][0],a[-1][1]],[1-1/n,1,1],fmt[not f],lw=lwd) # scruff
    else : 
        PLT.plot([a[-1][0],a[-1][1]],[1,1],fmt[not f],lw=lwd)
        PLT.plot([a[-1][0],a[-1][0]],[1,1-1/n],fmt[int(f)],lw=lwd)
    f = a[0].is_proper()
    if f : PLT.plot([a[0][0],a[0][1],a[0][1]],[0,0,1/n],fmt[int(f)],lw=lwd) # groin
    else : 
        PLT.plot([a[0][0],a[0][1]],[0,0],fmt[int(f)],lw=lwd)
        PLT.plot([a[0][1],a[0][1]],[0,1/n],fmt[not f],lw=lwd)
     
def Replot(a,fmt='w-',PLT=None) : 
    Plot(a,fmt=fmt,lwd=3,PLT=PLT)
    Plot(a,lwd=1,PLT=PLT)

def Dual(a):
    return([dual(q) for q in a])

def Sort(a) :
    L = [q[0] for q in a]
    L.sort()
    R = [q[1] for q in a]
    R.sort()
    return([k(L[i],R[i]) for i in range(len(a))])

def Shift(a,s) :
    return([s+A for A in a])

def Scale(a,s) :
    c = [s*A for A in a]
    return(Sort(c))

def Reciprocal(a) :
    return [1/A for A in a]
    
def Negate(a) :
    return(Sort([-A for A in a]))     # why is Sort necessary...it shouldn't be

def Condense(z,n=STEPS) :
    many = len(z)
    if many == n : return(z)
    sk = many / n
    return([z[i] for i in range(many) if i % sk == 0])

def Mix(a,b,wa=1,wb=1): # weighted stochastic mixture
    if wa < 1 or wb < 1 : 
        wa , wb = continued(wa/wb) # need whole numbers for the weights
    return Condense(Sort(a*wa+b*wb))

def Env(a,b) :
    return [A.env(B) for A,B in zip(a,b)]

def Imp(a,b) :
    return [A.imp(B) for A,B in zip(a,b)]

depdic = {'i':'independence', 'p':'perfect', 'o':'opposite', 'f':'Fréchet', '?':'Fréchet', '+': 'positive', '-': 'negative'}
depcol = {'i':'r', 'p':'g', 'o':'k', 'f':'b', '?':'b', '+': 'c', '-': 'o'}

def Add(a,b,d='f'):                          # default dependence should be 'f'
    if d=='i' : return(AddIndependent(a,b))
    if d=='p' : return(AddPerfect(a,b))
    if d=='o' : return(AddOpposite(a,b))
    if d=='f' : return(AddFrechet(a,b))
    #if d=='+' : return(AddPositive(a,b))
    #if d=='-' : return(AddNegative(a,b))
    
def AddIndependent(a,b) :      # sum under independence using Cartesian product
    c = []
    for A in a:
        for B in b:
            c = c + [A+B]
    return(Condense(Sort(c)))

def AddPerfect(a,b) :
    c = []
    for A,B in zip(a,b):
            c = c + [A+B]
    return(c)   # should be no need for Sort(c)

def AddOpposite(a,b) :
    c = []
    for A,B in zip(a,b[::-1]):  # [::-1] is the 'slicing' way to reverse a list
            c = c + [A+B]
    return(Sort(c))   # we need Sort(c)!

def AddFrechet(a,b) : 
    zu = [0.0] * STEPS;   xu = [q[0] for q in a];   yu = [q[0] for q in b]
    zd = [0.0] * STEPS;   xd = [q[1] for q in a];   yd = [q[1] for q in b]
    #s = ''; b = ' '
    for i in range(STEPS) :  # i in 1:STEPS
          #s = '';  
          outlier = float('inf') 
          for j in range(i,STEPS) :
            #s = s+str(j)+b+str(i-j+STEPS-1)+b
            here = xd[j] + yd[i-j+STEPS-1]
            if here<outlier : outlier = here
          #print(s)
          #s = ''; 
          zd[i] = outlier
          outlier = float('-inf')
          for j in range(i+1) :
            #s = s+str(j)+b+str(i-j)+b
            here = xu[j] + yu[i-j]
            if here>outlier : outlier = here
          #print(s)
          zu[i] = outlier   
    return([k(zu[i],zd[i]) for i in range(STEPS)])

def Subtract(a,b,d='f') :                    # default dependence should be 'f'
    if d=='i' : return(SubtractIndependent(a,b))
    if d=='p' : return(SubtractPerfect(a,b))
    if d=='o' : return(SubtractOpposite(a,b))
    if d=='f' : return(SubtractFrechet(a,b))
    #if d=='+' : return(SubtractPositive(a,b))
    #if d=='-' : return(SubtractNegative(a,b))
  
def SubtractIndependent(a,b) :          # differences using a Cartesian product
    c = []
    for A in a:
        for B in b:
            c = c + [A - B]
    return(Condense(Sort(c)))

def SubtractOpposite(a,b) :
    c = []
    for A,B in zip(a,b[::-1]):  # [::-1] is the 'slicing' way to reverse a list
        c = c + [A-B]
    return(c)   # should be no need Sort(c)

def SubtractPerfect(a,b) :     # Risk Calc actually gives the wrong answer here
    # see https://sites.google.com/site/riskcalcblog/corrections/subtraction
    c = []
    for A,B in zip(a,b):
        c = c + [A-B]
    return(Sort(c))
    # a = U([2,5], [6,7]); Plot(a)
    # c = Subtract(a,a,'p'); Plot(c)
    # d = Subtract(a,a,'o'); Plot(d)

def SubtractFrechet(a,b) :
    return(AddFrechet(a,Negate(b)))

#  if (op=='/') return(frechetconv.pbox(x,reciprocate.pbox(y),'*'))

# should it be able to just touch zero??

#  if (op=='*') if (straddlingzero(x) | straddlingzero(y)) return(imp(balchprod(x,y),naivefrechetconv.pbox(x,y,'*')))


def Multiply(a,b,d='f'):                          # default dependence should be 'f'
    if d=='i' : return(MultiplyIndependent(a,b))
    if d=='p' : return(MultiplyPerfect(a,b))
    if d=='o' : return(MultiplyOpposite(a,b))
    if d=='f' : return(MultiplyFrechet(a,b))
    #if d=='+' : return(MultiplyPositive(a,b))
    #if d=='-' : return(MultiplyNegative(a,b))
    
def MultiplyIndependent(a,b) :      # sum under independence using Cartesian product
    c = []
    for A in a:
        for B in b:
            c = c + [A * B]
    return(Condense(Sort(c)))

def MultiplyPerfect(a,b) :
    c = []
    for A,B in zip(a,b):
            c = c + [A * B]
    return(c)   # should be no need for Sort(c)

def MultiplyOpposite(a,b) :
    c = []
    for A,B in zip(a,b[::-1]):  # [::-1] is the 'slicing' way to reverse a list
            c = c + [A * B]
    return(Sort(c))   # we need Sort(c)!

def MultiplyFrechet(a,b) :    ###### this isn't working yet 
    
    #if either argument straddles zero, abort

    zu = [0.0] * STEPS;   xu = [q[0] for q in a];   yu = [q[0] for q in b]
    zd = [0.0] * STEPS;   xd = [q[1] for q in a];   yd = [q[1] for q in b]
    #s = ''; b = ' '
    for i in range(STEPS) :  # i in 1:STEPS
          outlier = float('inf') 
          for j in range(i,STEPS) :
            here = xd[j] * yd[i-j+STEPS-1]
            if here<outlier : outlier = here
          zd[i] = outlier
          outlier = float('-inf')
          for j in range(i+1) :
            here = xu[j] * yu[i-j]
            if here>outlier : outlier = here
          zu[i] = outlier   
    return([k(zu[i],zd[i]) for i in range(STEPS)])

def Red(a) : Plot(a,fmt='r-')
def Blue(a) : Plot(a,fmt='b-')
def Green(a) : Plot(a,fmt='g-')
def Purple(a) : Plot(a,fmt='m-')
def Cyan(a) : Plot(a,fmt='c-')
def Black(a) : Plot(a,fmt='k-')

def Plotabc(a,b,c,tt=None,PLT=None) :
    if PLT is None : PLT = plt
    Plot(a,lab='a',PLT=PLT); Plot(b,lab='b',PLT=PLT); 
    if not c is None : Plot(c,lab='c',PLT=PLT)
    if not tt is None : PLT.title(tt)

Uniform = U 
Exponential = E
Triangular = T
Interval = Minmax = I
plotbox = Plot

def Comp(a, b):
    c = eval(a)
    if not type(c)==Kaucher: c = Kaucher(c)
    ok = c == b
    if not ok : ok = k.dist(c,b) < Kaucher.negligibledifference
    if ok : print('{:>17}'.format(a),'  =  ',b)
    else : print(a,'\t',b, '\t',c, '******* ERROR *********')

###########################################################################
# # % % Generalized p-box constructors

# # Scott thinks these make sense
# Plot(U(k(1,2), k(3,5))); plt.show()        # ordinary p-box
# Plot(U(     1, k(3,5))); plt.show()        # ordinary p-box
# Plot(U(k(2,1), k(5,3))); plt.show()        # Kaucher p-box
# Plot(Dual(U(     1, k(3,5)))); plt.show()  # Kaucher p-box
# Plot(U(k(2,1), k(3,5))); plt.show()        # Kaucher p-box
# Plot(Dual(U(k(2,1), k(3,5)))); plt.show()  # Kaucher p-box

# Kaucher p-boxes are defined as stochastic mixtures of Kaucher intervals. 
# From such structures, Kaucher intervals can be randomly generated through 
# the inverse transform in the traditional way.  A random number from U(0,1) 
# specifies a pair of values x1 and x2 respectively from the left and right 
# edges of a p-box making a generalized interval [x1, x2].  Such intervals 
# from the first two example p-boxes 
#
#    Uniform(k(1,2), k(3,5))
#    Uniform( 1, k(3,5))
#
# are traditional proper intervals, where x1<=x2. The ones from the third p-box
#
#    Uniform([2,1], [5,3])) 
#    
# will result in improper intervals (i.e., x1>x2). The random intervals from 
# the fifth p-box 
#
#    Uniform( [2,1], [3,5])) 
#    
# will be intervals that are sometimes proper and sometimes improper.

###########################################################################
# % % Are these also generalized p-boxes?

# Plot(U(k(3,5), k(1,2))); plt.title('gibberish?'); plt.show()
# Plot(U(k(3,5),      1)); plt.title('gibberish?'); plt.show()
# Plot(U(k(5,3), k(2,1))); plt.title('gibberish?'); plt.show()
# Plot(U(k(3,5), k(2,1))); plt.title('gibberish?'); plt.show()

# Scott thinks these probably DON'T make sense; the uniform's minimum is larger 
# than its maximum.  Or could they make sense somehow in this bizarro Kaucher 
# world?

# Yan Wang notes: "At this point, I'm not sure about the significance of these 
# objects. My first feeling is that these might be useful purely arithmetically.  
# The most important aspect of a Kaucher p-box would be arithmetic convenience. 
# Improper intervals are "negative" intervals that significantly simplify the 
# arithmetic and computation, the same way negative numbers do with respect 
# to positive numbers. Negative numbers are abstract and were introduced 
# not long ago to make calculation easier. It would be great to keep open the 
# question about how Kaucher p-boxes can help reduce computational complexity 
# while designing them."

###########################################################################
# %% Generalized p-box constructions

def Subplot(where,a,fmt='',lab='',lwd='',tit=''):
    #n = len(a)
    if lwd == '' : lwd = 2
    if not tit == '' : where.set_title(tit)
    if not lab == '' : where.text(Mid(a),0.1,lab)
    Plot(a,fmt=fmt,lwd=lwd, PLT=where)

def someConstructions() :
    figure, axis = plt.subplots(3, 2)
    Subplot(axis[0,0],Uniform( [1,2],  [3,5] ),lab='U([1,2],[3,5])')
    Subplot(axis[0,1],Uniform(     1,  [3,5] ),lab='U(1,[3,5])')
    Subplot(axis[1,0],Uniform( [2,1],  [5,3] ),lab='U([2,1],[5,3])')
    Subplot(axis[1,1],Uniform(     1,  [5,3] ),lab='U(1,[5,3])')
    Subplot(axis[2,0],Uniform( [2,1],  [3,5] ),lab='U([2,1],[3,5])')
    Subplot(axis[2,1],Uniform( [1,2],  [5,3] ),lab='U([1,2],[5,3])')
    plt.show()

someConstructions()  

print("""
These are six UNIFORM Kaucher p-boxes constructed with U().

Kaucher p-boxes are defined as stochastic mixtures of Kaucher intervals. 
From such structures, Kaucher intervals can be randomly generated through 
the inverse transform in the traditional way.  A random number from U(0,1) 
specifies a pair of values x1 and x2 respectively from the left and right 
edges of a p-box making a generalized interval [x1, x2].  Such intervals 
from the first two example p-boxes 

   Uniform(k(1,2), k(3,5))
   Uniform( 1, k(3,5))

are traditional proper intervals, where x1<=x2. The ones from the third p-box

   Uniform([2,1], [5,3])) 
   
will result in improper intervals (i.e., x1>x2). The random intervals from 
the fifth p-box 

   Uniform( [2,1], [3,5])) 
   
will be intervals that are sometimes proper and sometimes improper.
""")

###########################################################################
# %% Gibberish constructions (are these also generalized p-boxes?)

def notsomeConstructions() :
    figure, axis = plt.subplots(2, 2)
    Subplot(axis[0,0],Uniform( [3,5], [1,2]),lab='U([3,5],[1,2])')
    Subplot(axis[0,1],Uniform( [3,5],     1),lab='U([3,5],1)')
    Subplot(axis[1,0],Uniform( [5,3], [2,1]),lab='U([5,3],[2,1])')
    Subplot(axis[1,1],Uniform( [3,5], [2,1]),lab='U([3,5],[2,1])')
    plt.show()
    
notsomeConstructions()  

print("""
Four 'GIBBERISH' abominations are condensed onto a multi-graph plot.  These
dubious constructions look like survival functions (i.e., complementary 
cumulative distributions or CCDFs) but they're not; they're just upside down.

Scott thinks these probably DON'T make sense; the uniform's minimum is larger 
than its maximum.  Or could they make sense somehow in this bizarro Kaucher 
world?

Yan Wang notes: "At this point, I'm not sure about the significance of these 
objects. My first feeling is that these might be useful purely arithmetically.  
The most important aspect of a Kaucher p-box would be arithmetic convenience. 
Improper intervals are "negative" intervals that significantly simplify the 
arithmetic and computation, the same way negative numbers do with respect 
to positive numbers. Negative numbers are abstract and were introduced 
not long ago to make calculation easier. It would be great to keep open the 
question about how Kaucher p-boxes can help reduce computational complexity 
while designing them."
""")

###########################################################################
# %% The normal distribution 

def P(w,pp) :
    col = ['c','k']
    Plot(N(*pp), PLT=w); 
    p = k(pp[0]); w.text(1.3,0.85,k.colorless(p),c=col[k.is_proper(p)])
    p = k(pp[1]); w.text(1.3,0.65,k.colorless(p),c=col[k.is_proper(p)])

figure, W = plt.subplots(4, 3, sharex=True, sharey=True)
P(W[0,0], (      8,       1))
P(W[1,0], ( k(8,9),       1))
P(W[2,0], (      8,  k(1,2)))
P(W[3,0], ( k(8,9),  k(1,2)))

P(W[0,1], ( k(8,9),  k(1,2)))
P(W[1,1], (~k(8,9),       1))
P(W[2,1], (      8,  ~k(1,2)))
P(W[3,1], (~k(5,9), ~k(1,2)))

P(W[0,2], ( k(8,9), ~k(1,2)))
P(W[1,2], (~k(5,9),  k(1,1)))
P(W[2,2], (~k(5,9), k(1,2)))
P(W[3,2], ( ~k(8,9),  ~k(1,2)))
plt.show()

print("""
Twelve generalized NORMAL p-boxes are shown.  The interval mean and interval 
standard deviation are shown in the upper, left corner of each graph.  The 
first column of graphs depict ordinary p-boxes defined by proper intervals, 
with red edges on the left.  Inverting the interval for the mean can yield 
an inverted p-box.  

Surprisily, merely inverting the interval for standard deviation doesn't 
seem to have any impact on the resulting p-box.  Compare, for instance, the 
second and third graphs in the top row, or the first and second on the third 
row, or the third graph on the third line with the second on the fourth. This
is due to the formulation used in the N() constructor.  However, it is by no
means clear that this is a reasonable one.

Leslie doesn't see anything wrong with this formulation of normal p-boxes.
He suggests that it seems reasonable that inverting the standard deviation
would have no effect on the p-box.
""")

###########################################################################
# %% More systematic exploration of Kaucher normal p-boxes 

mm = [k(12,13), k(10,13), k(12,13), k(10,13)]
ss = [k(1,2), k(1,2), k(1,3), k(1,3)]
figure, W = plt.subplots(4, 4, sharex=True, sharey=True)
for i in range(4) : P(W[i,0], (     mm[i],       ss[i]))
for i in range(4) : P(W[i,1], (dual(mm[i]),      ss[i]))
for i in range(4) : P(W[i,2], (     mm[i],  dual(ss[i])))
for i in range(4) : P(W[i,3], (dual(mm[i]), dual(ss[i])))
plt.show()

print("""
Again, the mean and standard deviation for each p-box are inscribed in the 
upper, left corner of each graph.  Reprising the observation that inverting 
the standard deviation does not seem to affect the p-box, the graphs in the 
first column are the same as those in the third, and those in the second are 
the same as those in the fourth.  As the observation calls into question the
formulation used in the N() constructor, we ask what checks or tests we might
employ to determine whether the formulation is correct.

If A = N(m,s), what are the parameters m' and s' such that N(m',s') is ~A?
The sum of independent normals should be normal; does this have an analog for 
Kaucher p-boxes?  And differences and sums with their own duals should give 
a result, at least under some dependence, that's comparable to the elementary 
self-combinations seen with intervals, e.g.,
        a = k(2,3)
        a - ~a      # k(0, 0)
        a + ~a      # k(5, 5)
A Kaucher p-box is a stochastic mixture of Kaucher intervals.

Leslie doesn't see anything wrong with this formulation of normal p-boxes.
He suggests that it seems reasonable that inverting the standard deviation
would have no effect on the p-box.
""")

###########################################################################
# %% The triangular distribution 

# The routine T() for constructing triangular p-boxes is not correct, and the
# earlier version T2() is worse.  T() gives proper p-boxes even with improper 
# parameters.  But it's not clear how the algorithm should behave.

def P(w,pp) :
    col = ['c','k']
    Plot(T(*pp),lwd=6, PLT=w); 
    Replot(T2(*pp), PLT=w) 
    p = k(pp[0]); w.text(1,0.85,k.colorless(p),c=col[k.is_proper(p)])
    p = k(pp[1]); w.text(1,0.65,k.colorless(p),c=col[k.is_proper(p)])
    p = k(pp[2]); w.text(1,0.45,k.colorless(p),c=col[k.is_proper(p)])

figure, W = plt.subplots(3, 3, sharex=True, sharey=True)
P(W[0,0], (      1,       3,      6))
P(W[0,1], (      1,  k(2,4),      6))
P(W[0,2], ( k(1,2),  k(3,4), k(5,6)))
P(W[1,0], (      1, ~k(2,4),      6))
P(W[1,1], (~k(1,2), ~k(3,4),      6))
P(W[1,2], (~k(1,3), ~k(2,4),      6))
P(W[2,0], (~k(1,3),  k(2,4),      6))
P(W[2,1], ( k(1,3), ~k(2,4),      6))
P(W[2,2], (~k(1,3), ~k(2,4), ~k(5,6)))
plt.show()

print("""
The formulation for TRIANGULAR p-box would seem to be even trickier.  Two
different formulations are proposed, and it may be that neither makes sense.
The results from T() are shown with thick lines.  The results from T2() are
shown with thin lines (on slightly thicker white background so the two sets 
of results are easily distinguishable).
""")

###########################################################################
# %% Directed rounding

def R(w,m,M) :
    a = U(m,M);         Plot(a,lwd=2,PLT=w);    
    A = U(m,M,10);      Plot(A,PLT=w)

figure, W = plt.subplots(3, 2, sharex=True, sharey=True)
R(W[0,0],1,4)       
R(W[1,0],k(1,2),k(3,4))
R(W[2,0],k(2,1),k(4,3))
R(W[0,1],k(1,2),k(4,3))
R(W[1,1],1,k(2,4))
R(W[2,1],1,k(4,2))
plt.show()

print("""
The theoretical distribution or p-box can be condensed to only ten levels.
Normally, OUTWARD-directed rounding is used to capture representation error.
But if the edges of the p-box are improper, I guess we need INWARD rounding.
""")

###########################################################################
# %% The basic interval intersection (imp) and envelope (env) operations

def iplot(x,  PLT=None) :
    if PLT is None : PLT = plt
    fmt = ['m-','-'][x.is_proper()]
    PLT.plot([x[0],x[0],x[1],x[1]], [0,1,1,0], fmt)

thr = 3
AA = [k(1,2), k(1,3), k(1,4)]
BB = [k(3,4), k(2,4), k(2,3)]
figure, W = plt.subplots(thr*2, thr*2, sharex=True, sharey=True)
for l in [0,1] :
    lt = l * thr    
    for j in [0,1] :
        jt = j * thr
        for i in range(thr) : 
            A = [AA[i], ~AA[i]][j % 2];        #a = Minmax(A)
            B = [BB[i], ~BB[i]][l % 2];        #b = Minmax(B)
            C = k.imp(A,B);                    #c = Imp(a,b)
            D = k.env(A,B);                    #d = Env(a,b)
            w = W[0+lt,i+jt]; iplot(A, PLT=w); #Plot(a,PLT=w)
            w = W[0+lt,i+jt]; iplot(B, PLT=w); #Plot(b,PLT=w)
            w  =W[1+lt,i+jt]; iplot(C, PLT=w); #Plot(c,PLT=w)
            w = W[2+lt,i+jt]; iplot(D, PLT=w); #Plot(d,PLT=w)
plt.show()

print("""
Several interval arithmetic INTERSECTION (imp) and ENVELOPE (env) operations.
Improper intervals are always shown in purple.  Inputs (A and B) are shown in 
the top row and in the fourth row.  The corresponding intersections imp(A,B) 
are show in the second and fifth rows.  The unions env(A,B) are shown below
them in the third and sixth rows.  If we use Minmax() constructors to make
Kaucher p-boxes analogous to these intervals, the intersections and envelopes 
exactly mirror the interval results.
""")  

###########################################################################
# %% Inverted p-boxes can arise from intersection calculations

def Q(w,wc,a,b) :
    Plot(a, PLT=w);                     Plot(b, PLT=w)
    Plot(a, PLT=wc, fmt='xkcd:gray');   Plot(b, PLT=wc, fmt='xkcd:gray')
    c = Imp(a,b);                       Plot(c, PLT=wc, lwd=2)

figure, W = plt.subplots(4, 4, sharex=True, sharey=True)
Q(W[0,0],W[1,0],N(k(5,6),1),N(k(9,10),1))
Q(W[0,1],W[1,1],N(k(5,6),1),N(~k(9,10),1))
Q(W[0,2],W[1,2],N(~k(5,6),1),N(k(9,10),1))
Q(W[0,3],W[1,3],N(~k(5,6),1),N(~k(9,10),1)) # intersection yields an "inverted envelope"
Q(W[2,0],W[3,0],N(k(5,10),1),N(k(8,9),1))
Q(W[2,1],W[3,1],U(k(1,3),k(10,12)),U(k(5,6),k(7,8)))
Q(W[2,2],W[3,2],U(~k(1,3),~k(10,12)),U(k(5,6),k(7,8))) # inverting second argument yields "inverted envelope"
Q(W[2,3],W[3,3],U(k(7,3),k(12,10)),U(k(3,5),k(6,9)))
plt.show()

print("""
Inverted p-boxes arise from INTERSECTING non-overlapping p-boxes.  Eight 
example intersections are shown.  In each graph of the top and third rows, 
there are two input p-boxes.  The corresponding intersection is shown in the 
graph below on the second and fourth rows.  These intersections are formed 
using the Imp() function according to the rule "from the first black to the 
second red".  That is, the intersection's black edge follows the leftmost 
black edge of a p-box, and its red edge follows the rightmost red edge of 
a p-box.  This is even true when the intersection is proper, as exemplified 
by the lower, left hand pair of graphs showing that, when the p-boxes do 
overlap, then the intersection is an ordinary (non-inverted) p-boxes.  
Intersections are from the first black to the second red, or in the multi-
argument case, to the last red.
""")

###########################################################################
# %% Enveloping

def Q(a,b,F,w=None,wc=None) :
    if w is None : w = plt
    if wc is None : wc = plt
    Plot(a, PLT=w);                     Plot(b, PLT=w)
    Plot(a, PLT=wc, fmt='xkcd:gray');   Plot(b, PLT=wc, fmt='xkcd:gray')
    c = F(a,b);                         Plot(c, PLT=wc, lwd=3)

figure, W = plt.subplots(4,4, sharex=True, sharey=True)
Q(  U(k(4,5),k(6,7)),   U(k(0,2),k(13,16)),  Env,W[0,0],W[1,0])
Q(  U(k(4,5),k(6,7)),   U(k(2,0),k(16,13)),  Env,W[0,1],W[1,1]) 
Q(  U(k(5,4),k(7,6)),   U(k(0,2),k(13,16)),  Env,W[0,2],W[1,2])
Q(  U(k(5,4),k(7,6)),   U(k(2,0),k(16,13)),  Env,W[0,3],W[1,3])
Q(  Minmax(k(5),k(1)),  Minmax(k(10),k(4)),  Env,W[2,0],W[3,0])
Q(  Minmax(k(5),k(1)),  Minmax(k(10),k(6)),  Env,W[2,1],W[3,1]) 
Q(  U(k(5,1),k(9,4)),   U(k(2,2),k(16,8)),  Env,W[2,2],W[3,2]) #*                 
Q(  U(k(5,1),k(9,4)),   U(k(3,0),k(16,13)),  Env,W[2,3],W[3,3]) #*
plt.show()

print("""
The mnemonic rule for ENVELOPING Kaucher p-boxes seems to be "first red to 
second black".  That is, the envelope's red edge is the leftmost red edge of 
the input p-boxes, and its black edge follows the rightmost black edge of the 
input p-boxes.  This is the opposite of the rule for Kaucher intersections.
Notice that the bounds of envelopes, like those of intersections, can be 
inverted, and the bounds may even cross one another.  We saw above that the 
intersection of inverted p-boxes can look like their envelope.  Conversely, 
we see among these examples that the envelope of inverted p-boxes can looks 
like their intersection (first graph of the last row), or what we might call 
their interior (last graph of the second row, and second of last).  The 
envelope of improper p-boxes may be either improper (e.g., last graph of 
second row), or proper (e.g., second graph of the last row).  Insofar as 
improper p-boxes don't overlap, the envelope is their interior space between 
them.  These behaviors extend what we see in envelopes of Kaucher intervals.
For instance, 
    env(k(4,1),k(5,3))  # k(4, 3)    # both improper and ranges don't overlap
    env(k(3,1),k(5,4))  # k(3, 4)    # both improper but ranges do overlap
""")

###########################################################################
# %% Shifting and scaling, including negation

def Threeplot(where,a,b,c,tit='',lab=['a','b','c',''], fmt='',lwd='',sup=False):
    def onebox(i,a,fmt='') :
        where.text((Mid(a)+left(left(a)))/2,0.25,lab[i])
        Plot(a, fmt=fmt, lwd=lwd, PLT=where)          
    if lwd == '' : lwd = 2
    if not tit == '' : where.set_title(tit) #,loc='left')
    if sup : 
        onebox(3,c,fmt)
        return
    onebox(0,a,fmt); onebox(1,b,fmt); onebox(2,c,fmt)

    
s1 = 6
s2 = -1
a = E([1,2]); b = U(k(-3,0),k(8,5))
da = Dual(a); db = Dual(b); 
figure, ax = plt.subplots(2, 2, sharex=True, sharey=True)
figure.suptitle('Shifts and scalings')
Threeplot(ax[0,0], a, Shift( a,s1), Scale( a,s2),'E([1,2])',       ['','+'+str(s1),str(s2)+'*']); plt.yticks([])
Threeplot(ax[1,0],da, Shift(da,s1), Scale(da,s2),'Dual',           ['','+'+str(s1),str(s2)+'*']); plt.yticks([])
Threeplot(ax[0,1], b, Shift( b,s1), Scale( b,s2),'U([-3,0],[8,5])',['','+'+str(s1),str(s2)+'*']); plt.yticks([])
Threeplot(ax[1,1],db, Shift(db,s1), Scale(db,s2),'Dual',           ['','+'+str(s1),str(s2)+'*']); plt.yticks([])
plt.show()

###########################################################################
# %% Additions

def someAdditions(a,b,dep,s=False,f='') :
    da = Dual(a); db = Dual(b); 
    figure, ax = plt.subplots(2, 2, sharex=True, sharey=True)
    if dep=='' : dep = list(depdic.keys())[0:4]
    else : figure.suptitle('Addition assuming '+depdic[dep])
    for d in dep :
        if s : f = depcol[d]
        Threeplot(ax[0,0], a, b,Add( a, b,d),'a + b = c',fmt=f,sup=s); plt.yticks([])
        Threeplot(ax[1,0], a,db,Add( a,db,d),'a + ~b = c',['a','~b ','c',''],f,sup=s); plt.yticks([])
        Threeplot(ax[0,1],da, b,Add(da, b,d),'~a + b = c',['~a ','b','c',''],f,sup=s); plt.yticks([])
        Threeplot(ax[1,1],da,db,Add(da,db,d),'~a + ~b = c',['~a ','~b ','c',''],f,sup=s); plt.yticks([])

a = Shift(E([1,2]),3) 
b = U([11,12], [13,15])
someAdditions(a,b,'p')
someAdditions(a,b,'o')
someAdditions(a,b,'i')
someAdditions(a,b,'f')
someAdditions(a,b,'',True)  # superimpose all dependencies
plt.show()

print("""
Here are FOUR QUADCHARTS, corresponding to four dependence assumptions, plus
a FIFTH QUADCHART summarizing the results from all the dependencies.  In each
figure, the upper, left plot shows addition of two ordinary p-boxes.  How does
the sum vary when the addends are replaced with their duals? The lower, right 
plot shows the addition of two Kaucher p-boxes that are duals of the addends.  
It gives exactly the same picture, but with the colors reversed.  So the sum 
of duals is the dual of the sum.  The other plots show additions of one 
traditional p-box with a Kaucher p-box.  Interestingly, if you switch which 
addend was replaced by its dual, you get the dual of the same result.  Could 
this possibly be a general result?
""")

###########################################################################
# %% Additions with a fatter 'a' box so the bounds don't cross themselves

a = U([11,13], [14,25])
b = Shift(E([1,2]),29)
someAdditions(a,b,'p')
someAdditions(a,b,'o')
someAdditions(a,b,'i')
someAdditions(a,b,'f')
someAdditions(a,b,'',True)  # superimpose all dependencies
plt.show()

print("""
These 4+1 quadcharts repeat the same addition calculations, but now we have 
used a fatter p-box for the addend a so the resulting p-boxes don't cross 
themselves as much.
""")

###########################################################################
# %% Subtraction is addition of the negated subtrahend

figure, ax = plt.subplots(2, 3, sharex=True, sharey=True)
figure.suptitle('Subtraction is addition of the negated subtrahend')

w = ax[0,0];           lab=['A ','B ','A - B']
a = U([4,7],[11,12]);  b = T([13,14],[15,16],[17,18]) ## SOMETHING WRONG WITH THE T CONSTRUCTOR!
c = Subtract(a,b);     d = Add(a, Scale(b,-1))  
Threeplot(w, a, b, c,'A-B', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',      fmt='y-',lwd=2, lab=lab); plt.yticks([])
# this is not much of a check, as the Frechet algorithm uses Add with Negation

w = ax[1,0];           lab=['','','  ~A-~B']
a = Dual(a);           b = Dual(b); 
c = Subtract(a,b);     d = Add(a, Negate(b))
Threeplot(w, a, b, c,'Dual(A) - Dual(B)', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',                  fmt='y-',lwd=2, lab=lab); plt.yticks([])

w = ax[0,1];           lab=['A ','B ','A /-/ B']
a = U([4,7],[11,12]);  b = T([13,14],[15,16],[17,18]) ## SOMETHING WRONG WITH THE T CONSTRUCTOR!
c = Subtract(a,b,'p'); d = Add(a, Scale(b,-1),'o')
Threeplot(w, a, b, c,'A/-/B', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',      fmt='y-',lwd=2, lab=lab); plt.yticks([])

w = ax[1,1];           lab=['','','  A/-/~B']
a = (a);               b = Dual(b); 
c = Subtract(a,b,'p'); d = Add(a, Negate(b),'o')
Threeplot(w, a, b, c,'A /-/ Dual(B)', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',                  fmt='y-',lwd=2, lab=lab); plt.yticks([])

# Independent
w = ax[0,2];           lab=['A ','D ','  A |-| D']
a = U([4,7],[11,12]);  b = Shift(Negate(E([.2,2])),-10)
c = SubtractIndependent(a,b); d = AddIndependent(a, Scale(b,-1))
Threeplot(w, a, b, c,'A|-|D', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',      fmt='y-',lwd=2, lab=['']*3); plt.yticks([])

w = ax[1,2];           lab=['','','~A|-|D']
a = Dual(a);           #b = Dual(b); 
c = SubtractIndependent(a,b); d = AddIndependent(a, Negate(b))
Threeplot(w, a, b, c,'Dual(A) |-| D', fmt='',  lwd=5, lab=lab); plt.yticks([])
Threeplot(w, a, b, d,'',              fmt='y-',lwd=2, lab=['']*3); plt.yticks([])

plt.show()

print("""
These are some preliminary checks of subtraction.  The gold lining on the 
plots illustrates checking that subtraction yields the same result as the 
addition with the negated subtrahend.  For both independence and Frechet,
the dependence assumption for the addition is the same as for the addition.
But, if the subtraction assumes perfect dependence, of course the addition 
must assume opposite dependence, and vice versa.
""") 

###########################################################################
# %% Subtractions 

def someSubtractions(a,b,dep,s=False,f='') :
    da = Dual(a); db = Dual(b); 
    figure, ax = plt.subplots(2, 2, sharex=True, sharey=True)
    if dep=='' : dep = list(depdic.keys())[0:4]
    else : figure.suptitle('Subtraction assuming '+depdic[dep])
    for d in dep :
        if s : f = depcol[d]
        Threeplot(ax[0,0], a, b,Subtract( a, b,d),'a - b = d',['a','b ','d',''],f,sup=s); plt.yticks([])
        Threeplot(ax[1,0], a,db,Subtract( a,db,d),'a - ~b = d',['a','~b ','d',''],f,sup=s); plt.yticks([])
        Threeplot(ax[0,1],da, b,Subtract(da, b,d),'~a - b = d',['~a ','b','d',''],f,sup=s); plt.yticks([])
        Threeplot(ax[1,1],da,db,Subtract(da,db,d),'~a - ~b = d',['~a ','~b ','d',''],f,sup=s); plt.yticks([])


a = U([11,12], [13,15])
b = Shift(Scale(E([0.2,2.5]),1/3),3) 
someSubtractions(a,b,'p')
someSubtractions(a,b,'o')
someSubtractions(a,b,'i')
someSubtractions(a,b,'f')
someSubtractions(a,b,'',True)  # superimpose all dependencies
plt.show()

print("""
These are 4+1 QUADCHARTS for subtraction under four dependence assumptions, 
with the 5th quadchart summarizing the results of all the dependencies.  In 
each, the upper, left plot shows subtraction of two ordinary p-boxes.  How 
does it vary when the inputs are replaced with their duals? The lower, right 
plot shows the difference of two Kaucher p-boxes that are duals of the inputs.  
It gives just the same picture as the upper, left corner, but with the colors 
reversed.  So the difference of duals is the dual of the difference. The other 
plots show subtractions with one traditional p-box and a Kaucher p-box. If you 
switch which input was replaced by its dual, the result is dual to the other 
result.  These are pretty clearly general patterns.  In the last quadchart 
comparing the results from the various dependencies, the Frechet result in the 
tightest in the dual-dual subtraction, which tracks with what we know about 
Frechet calculations in backcalculations. 
""")

###########################################################################
# %% Addition of a p-box with its dual 

def someDualAdds(w, a):  
    def show(a,fmt,lwd='') :
        Plot(a,fmt=fmt,lwd=lwd,PLT=w); 
    b = Dual(a)                      
    show(a,'')
    show(AddFrechet(a,b),fmt='b')
    show(AddIndependent(a,b),fmt='r')
    show(AddOpposite(a,b),fmt='k') #,lwd=3)
    show(AddPerfect(a,b),fmt='g',lwd=2)
    # show(Scale(AddOpposite(a,b),0.5),fmt='xkcd:gray')
    show(Scale(AddPerfect(a,b),0.5),fmt='xkcd:gray')

figure, w = plt.subplots(3, 3, sharex=True, sharey=True)
figure.suptitle('Addition of a p-box its dual')

someDualAdds(w[0,0], I(2,7))
someDualAdds(w[0,1], U([2,5],[4,7]))
someDualAdds(w[0,2], N(k(4,5),0.3))

someDualAdds(w[1,0], I(7,2))
someDualAdds(w[1,1], U([2,4],[3,7]))
someDualAdds(w[1,2], U([2,4],[6,5]))

someDualAdds(w[2,0], U([2,4],[5,4]))  # all dependencies the same??!
someDualAdds(w[2,1], Shift(Scale(E([0.2,0.9]),0.9),2))
someDualAdds(w[2,2], Dual(U([2,3],[7,4])))
plt.show()

print("""
Every plot shows the addition of a p-box with its dual.  This is something to 
explore because, for any interval A = [a1,a2], A + dual(A) = a1 + a2.  Is there
an analog of this for p-boxes?  In each plot, the input p-box is on the left.  
It is a proper p-box if its left edge is red and its right edge is black.  
Otherewise it is an improper p-box.  The sums of the p-box with its dual, under 
various dependence assumptions, are shown on the right in each plot.  The sum 
assuming perfect dependence is shown in green.  Under independence, in red.  
Frechet is shown in blue.  The sum under opposite dependence is shown in black.
The green distribution is shown as thick lines because we presume that a p-box
is perfectly dependent with its dual.  In the first column, the results under 
various dependence models all coincide.  Dividing the green perfect result by 
two creates the gray distribution superimposed on the input p-box.  Note that 
some examples yield a scalar sum (which appears as a vertical line), but 
usually the sum is precise distribution.  I don't know of any cases for which 
the dual-sum is a non-degenerate p-box.
""")

###########################################################################
# %% Difference between a p-box and its dual 

def someDualDiffs(w, a):  
    def show(a,fmt,lwd='') :
        Plot(a,fmt=fmt,lwd=lwd,PLT=w); 
    b = Dual(a)                      
    show(a,'')
    show(SubtractFrechet(a,b),fmt='b')
    show(SubtractFrechet(a,b),fmt='b') #,lwd=3)
    show(SubtractIndependent(a,b),fmt='r')
    show(SubtractOpposite(a,b),fmt='k') #,lwd=3)
    show(SubtractPerfect(a,b),fmt='g',lwd=2)

figure, w = plt.subplots(3, 3, sharex=True, sharey=True)
figure.suptitle('Difference between a p-box and its dual')

someDualDiffs(w[0,0], I(2,7))
someDualDiffs(w[0,1], U([2,5],[4,7]))
someDualDiffs(w[0,2], N(k(4,5),0.3))

someDualDiffs(w[1,0], I(7,2))
someDualDiffs(w[1,1], U([2,4],[4,7]))
someDualDiffs(w[1,2], U([2,4],[6,5]))

someDualDiffs(w[2,0], U([2,4],[5,4]))  # all dependencies the same??!
someDualDiffs(w[2,1], Shift(Scale(E([0.2,0.9]),0.9),2))
someDualDiffs(w[2,2], Dual(U([3,4],[8,6])))
plt.show()

print("""
Every plot shows the subtraction A - Dual(A).  If A is any interval, the 
difference A - dual(A) = 0.  What's the analog of this for p-boxes?  In each 
plot, the input p-box is now on the right.  If proper, red is the left edge,
and black on the right.  The differences between the p-box and its dual, under 
various dependence assumptions, are shown on the left in each plot.  In these 
examples, we've highlighted the green result which assumes perfect dependence 
between the p-box and its dual because it seems to always yield a scalar zero.
The difference under other dependence models are shown in other colors. In red 
for independence. Frechet is shown in blue.  The difference under opposite 
dependence is shown in black.  In the first column, the results under various 
dependence models all coincide.  
""")

###########################################################################
# %% Division of a p-box by its dual 

def someDualQuotients(w, a):  
    def show(a,fmt,lwd='') :
        Plot(a,fmt=fmt,lwd=lwd,PLT=w); 
    b = Dual(a)                      
    show(a,'')
    show(MultiplyFrechet(a,Reciprocal(b)),fmt='b')
    show(MultiplyFrechet(a,Reciprocal(b)),fmt='b') #,lwd=3)
    show(MultiplyIndependent(a,Reciprocal(b)),fmt='r')
    show(MultiplyOpposite(a,Reciprocal(b)),fmt='k') #,lwd=3)
    show(MultiplyPerfect(a,Reciprocal(b)),fmt='g',lwd=2)

figure, w = plt.subplots(3, 3, sharex=True, sharey=True)
figure.suptitle('Quotient of a p-box divided by its dual')

someDualQuotients(w[0,0], I(2,7))
someDualQuotients(w[0,1], U([2,5],[4,7]))
someDualQuotients(w[0,2], N(k(4,5),0.3))

someDualQuotients(w[1,0], I(7,2))
someDualQuotients(w[1,1], U([2,4],[4,7]))
someDualQuotients(w[1,2], U([2,4],[6,5]))

someDualQuotients(w[2,0], U([2,4],[5,4]))  # all dependencies the same??!
someDualQuotients(w[2,1], Shift(Scale(E([0.2,0.9]),0.9),2))
someDualQuotients(w[2,2], Dual(U([3,4],[8,6])))
plt.show()

print("""
We can compare these self-dual calculations with analogues with intervals. 
What is the dependence between something and its dual?  These calculations 
would seem to suggest it might be perfect, but opposite would seem to be a 
more reasonable assumption.  How do the results (shown as black distributions) 
from the previous calculations that assume opposite dependence.

For reference, the elementary rules for self-combinations for intervals are 
   Expression            yields
     A + A  		      2 * A
     A - A  		      contains 0
     A * A  		      A ** 2
     A / A  		      contains 1 if A does not contain 0
     A + dual(A) 	      sum(A)
     A - dual(A) 	      0
     A * dual(A) 	      prod(A), or 0 if A contains 0 !
     A / dual(A) 	      1\n
The analogous self-combinations with p-boxes are different.  For instance, 
subtracting the dual of a p-box from itself under yields zero but only under 
perfect dependence.
""")

###########################################################################
# %% 





print("""
SCOTT HAS NOT FINISHED EDITING THIS AND SUBSEQUENT CELLS

SUGGESTIONS FOR WHAT THEY SHOULD BE LIKE ARE MOST WELCOME
      """)




a = U([0,1],[1,2]); Plot(a)
z = SubtractFrechet(a, Dual(a)); Plot(z)
plt.plot([0,0],[0,1],':'); Plot(a,lwd=4); Plot(z,'b-',lwd=3)


a = E([.2,2]); Plot(a)
z = SubtractFrechet(a, Dual(a)); Plot(z)
plt.plot([0,0],[0,1],':'); Plot(a,lwd=4); Plot(z,'b-',lwd=3)


Plot(a,lwd=6); Plot(g)
f = AddFrechet(a,a)  # two quantities whose p-boxes have the same shape

g = BackAddFrechet(f, (a))

Plot(a,lwd=6); Plot(g)



print("""
SCOTT HAS NOT FINISHED EDITING THIS AND SUBSEQUENT CELLS

SUGGESTIONS FOR WHAT THEY SHOULD BE LIKE ARE MOST WELCOME
      """)





###########################################################################
# %% Backcaculations

# Perfect dependence

# Okay, it looks like it works perfectly when the dependence is perfect!
# Given a + b = c, we compute an estimate for b.  Call the estimate B.
# Then we put B back into the formula a+B whose sum we call ab.  The
# original c and the compute ab match perfectly.

# First, let's set up a, b  and their true sum
a = Shift(E([1,2]),3); Plot(a,lab='a')
b = U([11,12], [13,15]); Plot(b, lab='b')
c = Add(a,b, 'p'); Plot(c,lab='c')
Plotabc(a,b,c); plt.show()

# Now solve for the backcalculation for an unknown b
B = Subtract(c, Dual(a), 'p')

# Check whether the solution actually works
ab = Add(a, B, 'p')
Plot(B,'c',lwd=8); Plot(b,lab='b');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# The yellow solution ab coincides with the original c, and the backcalculated
# solution B shown in cyan coincides with the original b.  Success!

# Now solve the other backcalculation for an unknown a
A = Subtract(c, Dual(b), 'p')

# Check whether the solution actually works
ab = Add(A, b, 'p')
Plot(A,'c',lwd=8); Plot(a,lab='a');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# The yellow solution ab coincides with the original c, and the backcalculated
# solution A shown in cyan coincides with the original a.  Success again.

# Opposite dependence

# The backcalculation also seems to work pretty nicely when opposite 'o'
# is used as the dependence for the Add() functions and perfect 'p' 
# dependence is used for the Subtract() function!  So the formula must
# be something like 
#
#        if a \+\ b = c, then B = c /-/ ~a
#
a = Shift(E([1,2]),3); Plot(a,lab='a')
b = U([11,12], [13,15]); Plot(b, lab='b')
c = Add(a,b, 'o'); Plot(c,lab='c')   # add, opposite
Plotabc(a,b,c); plt.show()

# Now solve for the backcalculation.
B = Subtract(c, Dual(a), 'p')        # backcalc, PERFECT

# Check whether the solution actually works.
ab = Add(a, B, 'o')  # opposite
Plot(B,'c',lwd=8); Plot(b,lab='b');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# The yellow solution ab coincides with the original c, and the backcalculated
# solution B shown in cyan is close to agreeing with the original b.  Success!
#
# So, in this case, it seems that the constraints of c are indeed respected
# as required, but the backcalculated b does not find the original b exactly.
# This is not surprising, or bad, I think.

# But now solve the other backcalculation for an unknown a.
A = Subtract(c, Dual(b), 'p')  # backcalc, PERFECT

# Check whether the solution actually works
ab = Add(A, b, 'o')  # opposite
Plot(A,'c',lwd=8); Plot(a,lab='a');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# So, the yellow solution ab comes close but clearly MISSES the original c, and
# the backcalculated solution A shown in cyan likewise misses the original a.
# Just not sure if this is a success or not.  We know there is a true a such
# that, when multiplied by b under opposite dependence, will yield c exactly,
# so why can't this calculation find it?  Do I have a bug somewhere?  What am 
# missing?





# maybe it's something about the severity of the a p-box that's cause a problem


a = Add(I(-1,1),Shift(E([1,2]),3)); Plot(a,lab='a')
b = U([11,12], [13,15]); Plot(b, lab='b')
c = Add(a,b, 'o'); Plot(c,lab='c')   # add, opposite
Plotabc(a,b,c); plt.show()

# Now solve for the backcalculation.
B = Subtract(c, Dual(a), 'p')        # backcalc, PERFECT

# Check whether the solution actually works.
ab = Add(a, B, 'o')  # opposite
Plot(B,'c',lwd=8); Plot(b,lab='b');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# The yellow solution ab coincides with the original c, and the backcalculated
# solution B shown in cyan is close to agreeing with the original b.  Success!
#
# So, in this case, it seems that the constraints of c are indeed respected
# as required, but the backcalculated b does not find the original b exactly.
# This is not surprising, or bad, I think.

# But now solve the other backcalculation for an unknown a.
A = Subtract(c, Dual(b), 'p')  # backcalc, PERFECT

# Check whether the solution actually works
ab = Add(A, b, 'o')  # opposite
Plot(A,'c',lwd=8); Plot(a,lab='a');   Plot(ab,'y',lwd=8); Plot(c,lab='c'); 
plt.show()
# So, the yellow solution ab comes close but clearly MISSES the original c, and
# the backcalculated solution A shown in cyan likewise misses the original a.
# Just not sure if this is a success or not.  We know there is a true a such
# that, when multiplied by b under opposite dependence, will yield c exactly,
# so why can't this calculation find it?  Do I have a bug somewhere?  What am 
# missing?

























# subtracting a Kaucher p-box ... this should be backcalculation
a = U([11,12], [13,15]); Plot(a)
b = Dual(Shift(E([1,2]),3)); Plot(b)
c = Subtract(a,b); Plot(c)
Plotabc(a,b,c)

# if you dual the minuend instead
# all colors reverse
a = Dual(U([11,12], [13,15])); Plot(a)
b = Shift(E([1,2]),3); Plot(b)
c = Subtract(a,b); Plot(c)
Plotabc(a,b,c)

# now repeat with a fatter p-box for the variable a so the 
# p-box results don't cross themselves

# subtracting a Kaucher p-box ... this should be backcalculation
c = U([1,3], [4,15]); Plot(c)
a = Dual(Shift(E([1,2]),15)); Plot(a)
b = Subtract(c,a); Plot(b)
Plotabc(a,b,c)

b = Condense(b); Plot(b); plt.show()   # too many discretization levels
Plot(c)
z = Add(b,a)      # check the forward calculation
Plot(z,'-b')   # it's pretty close I guess
# are the discrepancies due to the dependence assumptions I'm ignoring?

# if you dual the minuend instead of the subtrahend, all the colors reverse
# but what does this operation represent???
c = Dual(U([1,3], [4,15])); Plot(c)
a = Shift(E([1,2]),15); Plot(a)
b = Subtract(c,a); Plot(b)
Plotabc(a,b,c)










###########################################################################
# %% Control solutions

# let's try the control solution, which is backcalcuate with skinny c and fat a
# we just swap what we called c with a and vice versa

a = U([1,3], [4,15]); Plot(a)
c = Dual(Shift(E([1,2]),15)); Plot(c)
# can't fit all the uncertainty of a into the tiny constraint of c
b = Subtract(c,a); Plot(b)
# so b is inverted...it's a control solution     OOPS!  NOT SO!
Plot(a); Plot(b); Plot(c); plt.show()
b = Condense(b); Plot(b); plt.show()   # too many discretization levels
Plot(c)
z = Add(b,a)      # check the forward calculation
Plot(z,'-b')   # does this seem correct?
# doesn't seem correct; it's not even inverted


"""
###########################################################################
# Self-combinations

# see SomeDualDiffs

def someSelfCombinations(a,w) :
    da = Dual(a);   
    bi = Subtract(a, da,d='i'); Red(bi)
    bp = Subtract(a, da,d='p'); Green(bp)
    bo = Subtract(a, da,d='o'); Black(bo)
    bf = Subtract(a, da,d='f'); Blue(bf)  
    figure, ax = plt.subplots(2, 2, sharex=True, sharey=True)
    figure.suptitle(' assuming '+depdic[dep])
    
a = U(k(11,12),k(13,15))
a = N(k(8,9),k(1,2))

# For intervals, the elementary rules for self-combinations are these:
#
# print('Kaucher interval self-combinations')
# print('Expression       yields')
# print(' A + A  		 2 * A')
# print(' A - A  		 contains 0')
# print(' A * A  		 A ** 2')
# print(' A / A  		 contains 1 if A does not contain 0')
# print(' A + dual(A) 	 sum(A)')
# print(' A - dual(A) 	 0')
# print(' A * dual(A) 	 prod(A), or 0 if A contains 0 !')
# print(' A / dual(A) 	 1\n')
#
# For instance, subtracting the dual of a p-box from itself under yields zero
# but only under perfect dependence.
#
"""

"""
def Plot2(a,fmt='',lab='',lwd='', PLT=None): # connect the dots WITHOUT considering representation error 
    if PLT is None : PLT = plt
    n = len(a)
    if lwd == '' : lwd = 2
    if not lab == '' : PLT.text(Mid(a),1.008,lab)
    if fmt == '' : fmt = ('r-','k-')
    if isinstance(fmt, str) : fmt = (fmt,fmt)   
    pp = [q/(n) for q in range(n)]   
    L = min(a[0])
    PLT.plot([L]+[q[0] for q in a]+[a[-1][1]],[0]+pp+[1],fmt[0],lw=lwd)    
    pp = [q/(n) for q in range(1,n+1)]
    PLT.plot([L,a[0][1]]+[q[1] for q in a],[0,0]+pp,fmt[1],lw=lwd)

"""


"""

frechetconv.pbox <- frechetconv.pbox <- function(x,y,op='+') {
  if (op=='-') return(frechetconv.pbox(x,negate.pbox(y),'+'))
  if (op=='/') return(frechetconv.pbox(x,reciprocate.pbox(y),'*'))
  if (op=='*') if (straddlingzero(x) | straddlingzero(y)) return(imp(balchprod(x,y),naivefrechetconv.pbox(x,y,'*')))
  x <- makepbox(x)
  y <- makepbox(y)
  zu <- zd <- rep(0,STEPS)
  for (i in 1:STEPS)
        {
          j <- i:STEPS
          k <- STEPS:i
          zd[[i]]  <- min(do.call(op,list(x@d[j],y@d[k])))  #zd[[i]] <- min(x@d[j] + y@d[k])
          j <- 1:i
          k <- i:1
          zu[[i]]  <- max(do.call(op,list(x@u[j],y@u[k])))  #zu[[i]] <- max(x@u[j] + y@u[k])
        }
  # mean
  ml <- -Inf
  mh <- Inf
  if (op %in% c('+','-')) {
    ml <- do.call(op, list(x@ml, y@ml))
    mh <- do.call(op, list(x@mh, y@mh))
    }
  # variance
  vl <- 0
  vh <- Inf
  if (op %in% c('+','-')) {
    #zv = env(xv+yv-2* sqrt(xv*yv), xv+yv+2* sqrt(xv*yv))
    # vh <- x@v+y@v - 2*sqrt(x@v*y@v)
    # vl <- x@v+y@v + 2*sqrt(x@v*y@v)
    }
  pbox(u=zu + 0, d=zd + 0, ml = ml, mh = mh, vl=vl, vh=vh, dids=paste(x@dids,y@dids))
  }


"""


"""
# incorrect discretization of quantiles

STEPS = 101 # number of discretization levels in a p-box

QQ = [q/(STEPS-1) for q in range(STEPS)]
bOt = 0.001
tOp = 0.999

def U(m:Kaucher, M:Kaucher, n=STEPS) : # note that STEPS is NOT evaluated, but the original value is used
    m = Kaucher(m)
    M = Kaucher(M)
    qq = [q/(n-1) for q in range(n)]
    return([k(m[0]+q*(M[0]-m[0]),m[1]+q*(M[1]-m[1])) for q in qq])

def E(m, n=STEPS):  # use E(1/lambda) for the other parameterization
    m = Kaucher(m)
    qq = [q/(n-1) for q in range(n)]
    e = [k(-m[0]*ln(1-q), -m[1]*ln(1-q)) for q in qq]
    e[-1] = k(max(e[-2][0],e[-2][1]))  # infinite tail hack, should use tOp and bOt
    return(e)

def N(m,s,n=STEPS) :
    Q = [q/(n-1) for q in range(n)]
    Q[0] = bOt
    Q[-1] = tOp
    L1 = sps.norm.ppf(Q,loc=left(m),scale=left(s))
    L2 = sps.norm.ppf(Q,loc=left(m),scale=right(s))
    R1 = sps.norm.ppf(Q,loc=right(m),scale=left(s))
    R2 = sps.norm.ppf(Q,loc=right(m),scale=right(s))
    return [k(min(L1[_],L2[_]), max(R1[_],R2[_])) for _ in range(n)]

def T2(m:Kaucher,mid:Kaucher,M:Kaucher, n=STEPS):
    def pr(p,pm,m,mid,M):
        if p<=pm : return m + (p * (M - m) * (mid - m))**0.5
        else :     return M - ((1.0 - p) * (M - m) * (M - mid))**0.5
    m = Kaucher(m)
    mid = Kaucher(mid)
    M = Kaucher(M)
    qq = [q/(n-1) for q in range(n)]
    pm = (mid-m) / (M-m)
    c = []
    for q in qq:
        c = c + [Kaucher(pr(q,pm[0],m[0],mid[0],M[0]),pr(q,pm[1],m[1],mid[1],M[1]))]        
    return(c)   # dividing pm[0] and pm[1] prolly doesn't make sense  

def T(m:Kaucher,mid:Kaucher,M:Kaucher, n=STEPS):
    def triangular(m : float, mid : float, M : float) :
        def pr(p,pm,m,mid,M):
            if p<=pm : return m + (p * (M - m) * (mid - m))**0.5
            else :     return M - ((1.0 - p) * (M - m) * (M - mid))**0.5
        pp = [q/(n-1) for q in range(n)]
        pm = (mid-m) / (M-m)
        c = []
        for q in qq:
            c = c + [Kaucher(pr(q,pm,m,mid,M),pr(q,pm,m,mid,M))]        
        return(c)
    m = Kaucher(m)
    mid = Kaucher(mid)
    M = Kaucher(M)
    # brute force strategy , ... but it's clearly WRONG for improper parameters
    tt = [Kaucher(math.inf,-math.inf)] * STEPS 
    for i in [0,1] :
        for j in [0,1] :
            for k in [0,1] :
                tt = Env(tt,triangular(m[i],mid[j],M[k]))
    return tt

def I(m:Kaucher, M=None, n=STEPS) : # this is basically the minmax function
    m = Kaucher(m)
    if M is None : M = m
    else : M = Kaucher(M)
    qq = [q/(n-1) for q in range(n)]
    return([k(m[0],M[1]) for q in qq])    # not sure this is what it should be!
    # Plot(I(k(1,3)))
    # Plot(I(k(10,3)))
    # Plot(I(k(1,4), k(10,3))) # Should this yield [1,3]?  It wouldn't seem so.
    # # But what should it be?
    # a = I(k(1,4), k(10,3))
    # b = I(k(1,4).intersection(k(10,3))); Plot(a,lwd=6); Replot(b)
    # b = I(k(1,4).env(k(10,3))); Plot(a,lwd=6); Replot(b)
        
def Plot2(a,fmt='',lab='',lwd='', PLT=None): # connect the dots WITHOUT outward-direction 
    if PLT is None : PLT = plt
    n = len(a)
    if lwd == '' : lwd = 2
    if not lab == '' : PLT.text(Mid(a),1.008,lab)
    pp = [q/(n-1) for q in range(n)]
    if fmt=='' : 
        PLT.plot([q[0] for q in a]+[a[-1][0],a[-1][1]],pp+[1,1],'-r',lw=lwd)
        PLT.plot([a[0][0],a[0][1]]+[q[1] for q in a],[0,0]+pp,'-k',lw=lwd)
    else:
        PLT.plot([q[0] for q in a]+[a[-1][0],a[-1][1]],pp+[1,1],fmt,lw=lwd)
        PLT.plot([a[0][0],a[0][1]]+[q[1] for q in a],[0,0]+pp,fmt,lw=lwd)
    
def Plot(a,fmt='',lab='',lwd='', PLT=None):
    if PLT is None : PLT = plt
    n = len(a)-1
    if fmt == '' : fmt = ('r-','k-')
    if isinstance(fmt, str) : fmt = (fmt,fmt)
    if lwd == '' : lwd = 1
    if not lab == '' : PLT.text(Mid(a),1.008,lab)
    #PLT.plot()
    for i in range(n) :
        x = [a[i][0],a[i][0],a[i+1][0]]
        p = [i/n,(i+1)/n,(i+1)/n]
        PLT.plot(x,p,fmt[0],lw=lwd)
    for i in range(n) :
        x = [a[i][1],a[i+1][1],a[i+1][1]]
        p = [i/n,i/n,(i+1)/n]
        PLT.plot(x,p,fmt[1],lw=lwd)
    f = a[-1].is_proper()
    PLT.plot([a[-1][0],a[-1][1]],[1,1],fmt[not f],lw=lwd) 
    f = a[0].is_proper()
    PLT.plot([a[0][0],a[0][1]],[0,0],fmt[f],lw=lwd)
     
def Replot(a,fmt='w-',PLT=None) : 
    Plot(a,fmt=fmt,lwd=3,PLT=PLT)
    Plot(a,lwd=1,PLT=PLT)
    
STEPS = 10  
b = N(5,1)
b = U(0,1,10)
plt.plot([0.5,0.5],[0,1]); plt.plot([0,1],[.5,.5]); Plot(b)

# the bull's eye is inside the 5th box; 
# it should be at the pinch point between the 5th and 6th boxes

"""