In [1]:
## Imports

import numpy as np
import sympy as sp
import scipy
from scipy import linalg
from scipy import integrate
from scipy import interpolate
import matplotlib.pyplot as plt
import matplotlib
import random
from numpy.random import seed
from matplotlib.patches import Rectangle
from functools import partial
import colorsys
import quadprog
import cvxopt 
import math
from timeit import default_timer as timer
from scipy.optimize import *
In [2]:
# s is the number of computed moments and there is one moment for every time-step
def distr2moms1D(CMESol,s):
    moms = np.zeros((s,CMESol.shape[1]))
    
    for t in range(CMESol.shape[1]):
        for i in range(CMESol.shape[0]):
            moms[0,t] += i*CMESol[i,t]
            
        for i in range(CMESol.shape[0]):
            for m in range(1,s):
                moms[m,t] += CMESol[i,t]*(i-moms[0,t])**(m+1)
                
    return moms

def moms2distr1D(moms):
    distr = np.zeros((moms.shape[0]+1,moms.shape[1]))
    
    A = np.zeros((moms.shape[0]+1, moms.shape[0]+1))
    b = np.ones(moms.shape[0]+1)
    
    for t in range(moms.shape[1]):
        
        for i in range(moms.shape[0]+1):

            A[0,i] = 1
            A[1,i] = i

            for j in range(2,moms.shape[0]+1):

                A[j,i] = (i-moms[0,t])**j

        for i in range(1,moms.shape[0]+1):
            b[i] = moms[i-1,t]
        #print(np.linalg.cond(A))
        distr[:,t] = np.linalg.solve(A,b)
        
    return A,distr

## Perform Quasi-Entropic Closure and Entropic Closure in a compareable manner
# Idea slightly different than for separate approaches: Choose Reconstruction Method and use equivalent Moment Extraction

def distrReconstruction(extendedMoms, numStates, closureMethod, optMethod):
    
    numMoments4Closure = extendedMoms.shape[0] - 1
    bounds = Bounds(np.zeros(numStates), np.ones(numStates))
    
    # Return all Equality constraints as a function
    def eqConstraintFunction(x, extendedMoms, numStates):
    
        numMoments4Closure = extendedMoms.shape[0] - 1
        # collect all equality constraints
        equalityConstraintsArray = np.zeros(1+numMoments4Closure)

        # norm criterion probSum -1 = 0
        equalityConstraintsArray[0] = np.sum(x) - extendedMoms[0]


        for i in range(1,1+numMoments4Closure):
            if (i == 1): # criterion to meet expected value
                equalityConstraintsArray[1] = np.sum(x*np.linspace(0,numStates-1,numStates)) - extendedMoms[1]
            else: # criterion for higher order centered moments (skewness, curtosis, ect...)
                equalityConstraintsArray[i] = np.sum(x*np.power((np.linspace(0,numStates-1,numStates) - extendedMoms[1]),i)) - extendedMoms[i]

        return equalityConstraintsArray

    eq_cons = {'type': 'eq',
               'fun' : partial(eqConstraintFunction, extendedMoms = extendedMoms, numStates = numStates)}

    ineq_cons = {'type': 'ineq',
                 'fun' : lambda x: x}
    
    if (szenarioName == "szenariosMany"):
        numIterations = 250
    elif (szenarioName == "szenariosFew"):
        numIterations = 25
    else:
        print("Invalid Szenario!")
    
    if (closureMethod == "Entropy"):            
        def optFunction(x):
            return np.sum(x*np.log(x))    
    elif (closureMethod == "QuasiEntropy"):
        def optFunction(x):
            n = x.shape[0]
            return np.sum(np.power((1/n-x),2))  
    else: print("Invalid Closure Type!")
    
    xInit = np.ones(numStates)/numStates
    optResult = minimize(optFunction, xInit, method=optMethod, constraints=[eq_cons, ineq_cons], bounds=bounds, options={'maxiter': numIterations})

    reconstructedDistr = optResult.x
    
    return reconstructedDistr 

def calcClosureMomentfromReconstructedDistr(reconstructedDistr, extendedMoms):
    
    numStates = reconstructedDistr.shape[0]
    
    matrixReconstructedDistr = np.zeros((numStates,1))
    matrixReconstructedDistr[:,0] = reconstructedDistr
    reconstructedMoms = distr2moms1D(matrixReconstructedDistr, extendedMoms.shape[0])
    
    # calculate the discrepancy between the moments of the reconstructed distribution and the objective moments given by extendedMoms
    
    # Version with absolute errors
    momentDiscrepancy = np.sum(np.power((extendedMoms[1:]-reconstructedMoms[:extendedMoms.shape[0]-1,0]),2))
    
    # Version with relative errors
    #momentDiscrepancy = np.sum(np.power((extendedMoms[1:]-reconstructedMoms[:extendedMoms.shape[0]-1,0])/extendedMoms[1:],2))

    return reconstructedMoms[-1], momentDiscrepancy
In [3]:
## Functions for the CME Solution

# Some propensities of unphysical states are included even though they are negative
# This should not be a problem since they can not be reached

def state2index(state, maxStateValues):
    index = 0
    
    for d in range(state.shape[0]):
        index += state[d]*np.prod(maxStateValues[0:d]+1)
        
    return int(index)

def index2state(index, maxStateValues, dim):
    tempIndex = index
    state = np.zeros(dim)
    
    for d in range(dim-1,-1,-1):
        state[d] = np.floor_divide(tempIndex, np.prod(maxStateValues[0:d]+1))
        tempIndex = np.mod(tempIndex, np.prod(maxStateValues[0:d]+1))
        
    return state

def constructCMEMatrix(V, propensitiesGetter, k, maxStateValues):
    CMEMatrix = np.zeros((np.prod(maxStateValues+1), np.prod(maxStateValues+1)))
    
    for index in range(np.prod(maxStateValues+1)):
        
        state = index2state(index, maxStateValues, V.shape[0])
        props = propensitiesGetter(state, k)

        
        for reac in range(V.shape[1]):
            destinationState = state + V[:,reac]
            
            if (np.all(destinationState <= maxStateValues) and np.all(destinationState >= 0)):
                destinationIndex = state2index(destinationState, maxStateValues)
                
                CMEMatrix[destinationIndex,index] += props[reac]
                CMEMatrix[index,index] -= props[reac]
                
    return CMEMatrix

def solveCMENum(initState, tEnd, V, propensitiesGetter, k, maxStateValues):
    
    CMEMatrix = constructCMEMatrix(V, propensitiesGetter, k, maxStateValues)
    
    pInitIndex = np.zeros(np.prod(maxStateValues+1))
    pInitIndex[state2index(initState, maxStateValues)] = 1.0
    
    pEndIndex = scipy.linalg.expm(CMEMatrix*tEnd).dot(pInitIndex)
    
    return pEndIndex

def solveSystem(initState, numStates, numTimeSteps, timeStepLength, V, propensitiesGetter, k, maxStateValues):
    CMESol = np.zeros((numStates,numTimeSteps))

    for tIndex in range(numTimeSteps):
        CMESol[:,tIndex] = solveCMENum(initState, timeStepLength*tIndex, V, propensitiesGetter, k, maxStateValues)
        
        
    def rhs(t,X,k):
        propensities = propensitiesGetter(X,k)
        result = np.zeros(X.shape)
        for r in range(V.shape[1]):
            result += V[:,r]*propensities[r]

            
        return result
        
    RRESol = scipy.integrate.solve_ivp(partial(rhs,k = k), (0,(numTimeSteps-1)*timeStepLength), initState, method='RK45')
    
    #plotCMEnRREResults(CMESol,RRESol.y, RRESol.t, numStates, numTimeSteps, timeStepLength, maxStateValues)
    
    return CMESol, RRESol
In [4]:
def momentRHSSystem(t, x, maxOrder, numStates, k, closureMethod, optMethod, systemID):
    
    m=np.zeros(maxOrder)
    
    global reconstructionDuration
    global reconstructionAccuracy
    global reconstructionTime
    global reconstructions

    extendedMoms = np.append(1.0, x) 

    ## Choose Moment Closure Method
    if (closureMethod == "Zero"):
        extendedMoms = np.append(extendedMoms, np.array([0.0]))
    else:
        start = timer()
        reconstructedDistr = distrReconstruction(extendedMoms, numStates, closureMethod, optMethod)
        end = timer()

        reconstructions = np.vstack((reconstructions,reconstructedDistr))

        closureMoment, momentDiscrepancy = calcClosureMomentfromReconstructedDistr(reconstructedDistr, extendedMoms)
        extendedMoms = np.append(extendedMoms, np.array([closureMoment])) 

        reconstructionDuration = np.append(reconstructionDuration, end-start)
        reconstructionAccuracy = np.append(reconstructionAccuracy, momentDiscrepancy)
        reconstructionTime = np.append(reconstructionTime, t)

        #print("Time: " + str(t))
        #print("Reconstruction Duration: " + str(end-start))
        #print("Reconstruction Accuracy: " + str(momentDiscrepancy))

    #Berechnung
    # System 1: 0 -> A; 2A -> 0
    if (systemID == 1):
        for e in range(1,maxOrder+1):
            summe=0.0
            for j in range (1, e+1):
                summe=summe+int(e-j!=1)*k[0]*scipy.special.comb(e, j)*extendedMoms[e-j]+math.pow(-1,j)*scipy.special.comb(e, j)*math.pow(2, j)*(int(e-j!=1)*k[1]*extendedMoms[1]*(extendedMoms[1]-1)*extendedMoms[e-j]+int(e-j!=0)*k[1]*(2*extendedMoms[1]-1)*extendedMoms[e-j+1]+k[1]*extendedMoms[e-j+2])
            summe=summe+e*int(e>=3)*extendedMoms[e-1]*(-k[0]+2*(k[1]*extendedMoms[1]*(extendedMoms[1]-1)+k[0]*extendedMoms[2]))
            m[e-1]=summe
        return m

    # System 2: A+B -> 2A; A+B -> 2B
    elif (systemID == 2):
        for e in range(1,maxOrder+1):
            summe=0.0
            for j in range (1, e+1):
                summe=summe+(k[1]+(-1)**j*k[0])*scipy.special.comb(e, j)*(int(e-j!=1)*extendedMoms[1]*((numStates-1)-extendedMoms[1])*extendedMoms[e-j]+int(e-j!=0)*(-2*extendedMoms[1]+numStates-1)*extendedMoms[e-j+1]-extendedMoms[e-j+2])
            summe=summe+(k[0]-k[1])*e*int(e>=3)*extendedMoms[e-1]*(extendedMoms[1]*(numStates-1-extendedMoms[1])-extendedMoms[2])
            m[e-1]=summe

    # System 3: 2A -> B; B -> 2A
    elif (systemID == 3):
        for e in range(1, maxOrder+1):
            summe = 0.0
            for j in range(1, e+1):
                summe = summe + scipy.special.comb(e, j)*2**j*((-1)**j*(int(e-j!=1)*k[0]*extendedMoms[1]*(extendedMoms[1]-1)*extendedMoms[e-j] + int(e-j!=0)*k[0]*(2*extendedMoms[1]-1)*extendedMoms[e-j+1] + int(e-j!=-1)*k[0]*extendedMoms[e-j+2]) + int(e-j!=1)*0.5*k[1]*(numStates-1-extendedMoms[1])*extendedMoms[e-j] - 0.5*int(e-j!=0)*k[1]*extendedMoms[e-j+1])

            summe = summe + 2*e*int(e>=3)*extendedMoms[e-1]*(k[0]*(extendedMoms[1]*(extendedMoms[1]-1) + extendedMoms[2]) - 0.5*k[1]*(numStates-1-extendedMoms[1]))
            m[e-1]=summe

    return m
In [5]:
def runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,solMethod, optMethod):
        
    global reconstructionDuration
    global reconstructionAccuracy
    global reconstructionTime
    global reconstructions

    reconstructionDuration = np.array([0])
    reconstructionAccuracy = np.array([0])
    reconstructionTime = np.array([0])
    reconstructions = np.zeros((1,numStates))
    
    maxStateValues = np.array([numStates-1])
    
    initMoms = np.zeros(numMoments4Closure)
    initMoms[0] = initState

    # System 1: 0 -> A; 2A -> 0
    if (systemID == 1):
        V = np.array([[1,-2]])
        def propensitiesGetter(X,k):
            return np.array([k[0], 
                             X[0]*(X[0]-1)*k[1]])
        
    # System 2: A+B -> 2B; A+B -> 2A
    elif (systemID == 2):
        V = np.array([[-1,1]])
        def propensitiesGetter(X,k):
            return np.array([X[0]*(maxStateValues[0]-X[0])*k[0], 
                             X[0]*(maxStateValues[0]-X[0])*k[1]])
        
    # System 3: 2A -> B; B -> 2A
    elif (systemID == 3):
        V = np.array([[-2,2]])
        def propensitiesGetter(X,k):
            return np.array([X[0]*(X[0]-1)*k[0], 
                            (maxStateValues[0]-X[0])/2*k[1]])
        
        
    if (solMethod == "CME"):
        sol,_ = solveSystem(initState, numStates, numTimeSteps, timeStepLength, V, propensitiesGetter, k, maxStateValues)

    elif (solMethod == "QuasiEntropy"):
        arguments = (numMoments4Closure, numStates, k, "QuasiEntropy", optMethod, systemID)
        sol = scipy.integrate.solve_ivp(lambda t, x: momentRHSSystem(t, x, *arguments), (0,(numTimeSteps-1)*timeStepLength), initMoms, method="RK45",)
    
    elif (solMethod == "Entropy"):
        arguments = (numMoments4Closure, numStates, k, "Entropy", optMethod, systemID)
        sol = scipy.integrate.solve_ivp(lambda t, x: momentRHSSystem(t, x, *arguments), (0,(numTimeSteps-1)*timeStepLength), initMoms, method="RK45",)
    
    elif (solMethod == "Zero"):
        arguments = (numMoments4Closure, numStates, k, "Zero", optMethod, systemID)
        sol = scipy.integrate.solve_ivp(lambda t, x: momentRHSSystem(t, x, *arguments), (0,(numTimeSteps-1)*timeStepLength), initMoms, method="RK45")
    
    else:
        print("Invalid Solution Method")
    
    return sol, reconstructionDuration, reconstructionAccuracy, reconstructionTime, np.transpose(reconstructions)
In [6]:
def visInequalities(mean, var, maxValue):
    
    f = plt.figure()
    f.set_figwidth(10)
    f.set_figheight(5)
    
    x = np.linspace(0,maxValue,1000)
    xleft = np.linspace(0,mean,1000)
    xright = np.linspace(mean,maxValue,1000)

    plt.plot(x,1-mean/x, color = [0.8,0,0],linewidth = 3.0,  label = "Markov")
    plt.plot(x, (maxValue-mean)/(maxValue-x), color = [0.8,0,0],linewidth = 3.0)
    plt.plot(xright,1-var/np.power((xright-mean),2), color = [0,0,0.8],linewidth = 3.0,  label = "Chebychev")
    plt.plot(xleft,var/np.power(mean-xleft,2), color = [0,0,0.8],linewidth = 3.0)
    
    plt.fill_between(x, 1-mean/x, np.zeros(1000), color=[0.8,0,0], alpha=0.3)
    plt.fill_between(x, (maxValue-mean)/(maxValue-x), np.ones(1000), color=[0.8,0,0], alpha=0.3)
    plt.fill_between(xright, 1-var/np.power((xright-mean),2), np.zeros(1000), color=[0,0,0.8], alpha=0.3)
    plt.fill_between(xleft, var/np.power(mean-xleft,2), np.ones(1000), color=[0,0,0.8], alpha=0.3)
    
    plt.title(r"Bounds for $F_X$")
    plt.axis([0,maxValue,0,1])
    plt.legend()
    plt.show()
    

# check if there exists a distribution on the support [0,numStates-1] with the provided moments
# returns a vector with one entry for each moment inequality
def checkMomentFeasibility(moments, numStates):
    
    c = numStates-1
    E = moments[0]    
    
    feas = np.ones(moments.shape[0])
    
    if ((E < 0) or (E > c)):
        feas[0] = 0
    
    if (moments.shape[0] > 1):
        V = moments[1]
        if ((V > c*E - E**2) or (V < 0)):
            feas[1] = 0
            
    if (moments.shape[0] > 2):
        S = moments[2]
        if ((S > (c-3*E)*(V+E**2) + 2*E**3) or (S < -(c-3*E)*(V+E**2) - 2*E**3)):
            feas[2] = 0
    
    if (moments.shape[0] > 3):
        C = moments[3]
        if ((C > (c-4*E)*(S-2*E**3) + ((c-4*E)*3*E + 6*E**2)*(V+E**2) - 3*E**4) or (C < 0)):
            feas[3] = 0
            
    return feas
In [13]:
def plotMomentTimelines(solEy, solQy, solZy, solEt, solQt, solZt, momsC):
    
    %matplotlib inline
    matplotlib.rcParams.update({'font.size': 13})

    colormapE = np.array([colorEInFeas, colorE])
    colormapQ = np.array([colorQInFeas, colorQ])

    feasArrayE = np.zeros((solEy.shape[0], solEy.shape[1]), dtype = "int")
    feasArrayQ = np.zeros((solQy.shape[0], solQy.shape[1]), dtype = "int")

    for i in range(solEy.shape[1]):
        feasArrayE[:,i] = checkMomentFeasibility(solEy[:,i], numStates)

    for i in range(solQy.shape[1]):
        feasArrayQ[:,i] = checkMomentFeasibility(solQy[:,i], numStates)

    momentNames = [r"Mean $\mathbb{E}$", r"Variance $\mathbb{V}$", r"Skewness $\mathbb{S}$", r"Curtosis $\mathbb{C}$"]
    numXTicks = [6,5,3,2]

    for i in range(momsC.shape[0]):
        f = plt.figure()
        f.set_figwidth(10-3.4*i)
        f.set_figheight(5-1.7*i)

        plt.plot(np.linspace(0,timeStepLength*numTimeSteps,numTimeSteps),momsC[i,:], color = colorC,linewidth = 3.0, label = "CME")

        #plt.plot(solZt, solZy[i,:], color = colorZ,linewidth = 3.0, label = "Zero")
        plt.plot(solEt, solEy[i,:], color = colorE,linewidth = 1.0, label = "Entropy")
        plt.plot(solQt, solQy[i,:], color = colorQ,linewidth = 1.0, label = "Quasi-Entropy")
        
        plt.scatter(solEt, solEy[i,:], color = colormapE[feasArrayE[i,:]], linewidth = 1.0)
        plt.scatter(solQt, solQy[i,:], color = colormapQ[feasArrayQ[i,:]], linewidth = 1.0)
        
        plt.xticks(np.linspace(0,1,numXTicks[i]))
        #plt.title(momentNames[i],loc='left')
        
        #if (i == 0):
            #plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left',ncol=4)    
            #plt.xlabel(r"time $t$")
        #else:
            #plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        
        plt.axis([-0.05*(numTimeSteps*timeStepLength), 1.05*(numTimeSteps*timeStepLength), 0, np.amax(momsC[i,:]) + 0.4*(np.amax(momsC[i,:])-np.amin(momsC[i,:]))])
        #plt.axis([-0.05*(numTimeSteps*timeStepLength), 1.05*(numTimeSteps*timeStepLength), np.amin(momsC[i,:])-0.2*(np.amax(momsC[i,:])-np.amin(momsC[i,:])), np.amax(momsC[i,:]) + 0.4*(np.amax(momsC[i,:])-np.amin(momsC[i,:]))])
        plt.show()
        
def plotReconstructionStats(reconstructionTimeQ, reconstructionTimeE, reconstructionAccuracyQ, reconstructionAccuracyE, reconstructionDurationQ, reconstructionDurationE):
    %matplotlib inline
    matplotlib.rcParams.update({'font.size': 13})
    # delete the first netry that is 0 anyway
    reconstructionTimeQ = reconstructionTimeQ[1:]
    reconstructionTimeE = reconstructionTimeE[1:] 
    reconstructionAccuracyQ = reconstructionAccuracyQ[1:] 
    reconstructionAccuracyE = reconstructionAccuracyE[1:]
    reconstructionDurationQ = reconstructionDurationQ[1:] 
    reconstructionDurationE = reconstructionDurationE[1:]
    
    f = plt.figure()
    f.set_figwidth(10)
    f.set_figheight(4)

    plt.subplot(2, 1, 1)

    plt.scatter(reconstructionTimeQ, np.log10(reconstructionAccuracyQ), color = colorQ, s = 20, label = "Quasi-Entropy")
    plt.scatter(reconstructionTimeE, np.log10(reconstructionAccuracyE), color = colorE, s = 20, label = "Entropy")

    plt.title("Reconstruction Accuracy")
    #plt.axis([0,timeStepLength*numTimeSteps,0,numStates])
    plt.legend()

    plt.subplot(2, 1, 2)

    plt.scatter(reconstructionTimeQ, np.log10(reconstructionDurationQ), color = colorQ, s = 20, label = "Quasi-Entropy")
    plt.scatter(reconstructionTimeE, np.log10(reconstructionDurationE), color = colorE, s = 20, label = "Entropy")

    plt.title("Reconstruction Duration")
    #plt.axis([0,timeStepLength*numTimeSteps,0,numStates])
    plt.legend()
    plt.show()
    
    
    logMinReconstrDur = min(np.amin(np.log10(reconstructionDurationQ)), np.amin(np.log10(reconstructionDurationE)))
    logMaxReconstrDur = max(np.amax(np.log10(reconstructionDurationQ)), np.amax(np.log10(reconstructionDurationE)))
    
    f = plt.figure()
    f.set_figwidth(6)
    f.set_figheight(3)
    
    plt.hist(np.log10(reconstructionDurationQ), bins = np.linspace(logMinReconstrDur-np.power(10.0,-10), logMaxReconstrDur+np.power(10.0,-10),60), color = colorQ, label = "Quasi-Entropy")
    plt.hist(np.log10(reconstructionDurationE), bins = np.linspace(logMinReconstrDur-np.power(10.0,-10), logMaxReconstrDur+np.power(10.0,-10),60), color = colorE, label = "Entropy")
    plt.xticks([-2,-1.5,-1,-0.5,0], [r"$10^{-2}$", r"$10^{-1.5}$", r"$10^{-1}$", r"$10^{-0.5}$", r"$10^{0}$"])
    
    #plt.title(r"Duration of Single Distr. Reconstruction")
    #plt.ylabel(r"frequency $[1]$")
    #plt.xlabel(r"duration $[s]$")
    #plt.legend(loc = "upper center")
    plt.show()
    
def plotReconstructions(reconstructionTimeE, reconstructionTimeQ, reconstructionTimeC, reconstructionsE, reconstructionsQ, reconstructionsC):
    %matplotlib inline
    matplotlib.rcParams.update({'font.size': 13})
    f = plt.figure()
    f.set_figwidth(20)
    f.set_figheight(4)
    #f.suptitle("Reconstructions for Entropy and Quasi-Entropy Closure vs. CME Solution")
    
    plt.subplot(1, 3, 1)
    
    coarseGrid = np.linspace(0,reconstructionsE.shape[0]-1,reconstructionsE.shape[0])
    fineGrid = np.linspace(0,reconstructionsE.shape[0]-1,reconstructionsE.shape[0]*100)

    for i in range(reconstructionsE.shape[1]):
        spline = scipy.interpolate.PchipInterpolator(coarseGrid,reconstructionsE[:,i])      
        plt.plot(fineGrid,spline(fineGrid), color=(1-reconstructionTimeE[i])*np.array(colorE),linewidth = 0.5, alpha = 1/10)
        plt.fill_between(fineGrid, spline(fineGrid),0, color=(1-reconstructionTimeE[i])*np.array(colorE), alpha= 1/20)
    
    #plt.ylabel(r"Probability $P(x)$")
    #plt.xlabel(r"States $x$")
    
    plt.subplot(1, 3, 2)
    
    coarseGrid = np.linspace(0,reconstructionsQ.shape[0]-1,reconstructionsQ.shape[0])
    fineGrid = np.linspace(0,reconstructionsQ.shape[0]-1,reconstructionsQ.shape[0]*100)

    for i in range(reconstructionsQ.shape[1]):
        spline = scipy.interpolate.PchipInterpolator(coarseGrid,reconstructionsQ[:,i])      
        plt.plot(fineGrid,spline(fineGrid), color=(1-reconstructionTimeQ[i])*np.array(colorQ),linewidth = 0.5, alpha = 1/10)
        plt.fill_between(fineGrid, spline(fineGrid),0, color=(1-reconstructionTimeQ[i])*np.array(colorQ), alpha= 1/20)
    
    #plt.xlabel(r"States $x$")
    
    
    plt.subplot(1, 3, 3)
    
    coarseGrid = np.linspace(0,reconstructionsC.shape[0]-1,reconstructionsC.shape[0])
    fineGrid = np.linspace(0,reconstructionsC.shape[0]-1,reconstructionsC.shape[0]*100)

    for i in range(reconstructionsC.shape[1]):
        spline = scipy.interpolate.PchipInterpolator(coarseGrid,reconstructionsC[:,i])      
        plt.plot(fineGrid,spline(fineGrid), color = (1-reconstructionTimeC[i])*np.array(colorC),linewidth = 0.5, alpha = 1/2.5)
        plt.fill_between(fineGrid, spline(fineGrid),0, color=(1-reconstructionTimeC[i])*np.array(colorC), alpha= 1/5)    
        
    #plt.xlabel(r"States $x$")
    
    plt.show()
In [8]:
def parameterReturner(szenarioNumber):
    
    if (szenarioNumber == 3):
        systemID = 1
        numMoments4Closure = 3
        k = np.array([1.0,1.0])
        return systemID, numMoments4Closure, k
    
    elif (szenarioNumber == 4):
        systemID = 2
        numMoments4Closure = 2
        k = np.array([2.25,1.0])
        return systemID, numMoments4Closure, k
    
    elif (szenarioNumber == 5):
        systemID = 2
        numMoments4Closure = 3
        k = np.array([2.25,1.0])
        return systemID, numMoments4Closure, k
    
    elif (szenarioNumber == 9):
        systemID = 3
        numMoments4Closure = 3
        k = np.array([1.0,2.0])
        return systemID, numMoments4Closure, k
    
    else:
        print("No Valid Szenario!")
        return 0
In [9]:
def simulateSystem4DiffClosures(szenarioNumber, szenarioName):
    
    print("Simulation " + str(szenarioNumber))
    
    systemID, numMoments4Closure, k = parameterReturner(szenarioNumber)
    
    solZ, reconstructionDurationZ, reconstructionAccuracyZ, reconstructionTimeZ, reconstructionsZ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"Zero", "SLSQP")
    print("Zero finished")
    
    solQ, reconstructionDurationQ, reconstructionAccuracyQ, reconstructionTimeQ, reconstructionsQ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"QuasiEntropy", "SLSQP")
    print("QuasiEntropy finished")

    solE, reconstructionDurationE, reconstructionAccuracyE, reconstructionTimeE, reconstructionsE = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"Entropy", "SLSQP")  
    print("Entropy finished")

    solC, _,_,_,_ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"CME", "SLSQP")
    print("CME finished")
    
    momsC = distr2moms1D(solC,numMoments4Closure)
    
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solZy" + str(szenarioNumber) + ".csv", solZ.y, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solQy" + str(szenarioNumber) + ".csv", solQ.y, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solEy" + str(szenarioNumber) + ".csv", solE.y, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solC" + str(szenarioNumber) + ".csv", solC, delimiter = ",")
    
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solZt" + str(szenarioNumber) + ".csv", solZ.t, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solQt" + str(szenarioNumber) + ".csv", solQ.t, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/solEt" + str(szenarioNumber) + ".csv", solE.t, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/momsC" + str(szenarioNumber) + ".csv", momsC, delimiter = ",")
    
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionsQ" + str(szenarioNumber) + ".csv", reconstructionsQ, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionsE" + str(szenarioNumber) + ".csv", reconstructionsE, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionTimeQ" + str(szenarioNumber) + ".csv", reconstructionTimeQ, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionTimeE" + str(szenarioNumber) + ".csv", reconstructionTimeE, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionAccuracyQ" + str(szenarioNumber) + ".csv", reconstructionAccuracyQ, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionAccuracyE" + str(szenarioNumber) + ".csv", reconstructionAccuracyE, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionDurationQ" + str(szenarioNumber) + ".csv", reconstructionDurationQ, delimiter = ",")
    np.savetxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionDurationE" + str(szenarioNumber) + ".csv", reconstructionDurationE, delimiter = ",")
In [10]:
def plotSystem4DiffClosures(szenarioNumber, szenarioName):
    
    solZy = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solZy" + str(szenarioNumber) + ".csv", delimiter = ",")
    solQy = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solQy" + str(szenarioNumber) + ".csv", delimiter = ",")
    solEy = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solEy" + str(szenarioNumber) + ".csv", delimiter = ",")
    solC = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solC" + str(szenarioNumber) + ".csv", delimiter = ",")
    
    solZt = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solZt" + str(szenarioNumber) + ".csv", delimiter = ",")
    solQt = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solQt" + str(szenarioNumber) + ".csv", delimiter = ",")
    solEt = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/solEt" + str(szenarioNumber) + ".csv", delimiter = ",")
    momsC = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/momsC" + str(szenarioNumber) + ".csv", delimiter = ",")
    
    reconstructionsQ = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionsQ" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionsE = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionsE" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionTimeQ = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionTimeQ" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionTimeE = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionTimeE" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionAccuracyQ = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionAccuracyQ" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionAccuracyE = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionAccuracyE" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionDurationQ = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionDurationQ" + str(szenarioNumber) + ".csv", delimiter = ",")
    reconstructionDurationE = np.loadtxt(szenarioName + "/" + str(szenarioNumber) + "/reconstructionDurationE" + str(szenarioNumber) + ".csv", delimiter = ",")

    plotMomentTimelines(solEy, solQy, solZy, solEt, solQt, solZt, momsC)
    plotReconstructionStats(reconstructionTimeQ, reconstructionTimeE, reconstructionAccuracyQ, reconstructionAccuracyE, reconstructionDurationQ, reconstructionDurationE)
    plotReconstructions(reconstructionTimeE, reconstructionTimeQ, np.linspace(0,timeStepLength*numTimeSteps,numTimeSteps), reconstructionsE, reconstructionsQ, solC)
In [11]:
timeStepLength = 0.01
initState = np.array([10])
numStates = 21
numTimeSteps = 100

colorC = (255/255, 192/255, 0/255,1.0)
colorZ = (80/255, 80/255, 80/255,1.0)
colorE = (0/255, 112/255, 192/255,1.0)
colorEInFeas = (0.0,0.0,0.0,1.0)
colorQ = (192/255, 0/255, 0/255,1.0)
colorQInFeas = (0.0,0.0,0.0,1.0)
In [14]:
szenarioNames = ["szenariosFew"]

for szenarioName in szenarioNames:
    for i in [3]:
        #simulateSystem4DiffClosures(i, szenarioName)
        plotSystem4DiffClosures(i, szenarioName)
In [ ]:
## Compare system runs with different closures
systemID = 1

timeStepLength = 0.01
initState = np.array([10])
numStates = 21
numMoments4Closure = 4
k = np.array([100.0,1.0])
numTimeSteps = 100

#solZ, reconstructionDurationZ, reconstructionAccuracyZ, reconstructionTimeZ, reconstructionsZ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"Zero", "SLSQP")
#solQ, reconstructionDurationQ, reconstructionAccuracyQ, reconstructionTimeQ, reconstructionsQ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"QuasiEntropy", "SLSQP")
#solE, reconstructionDurationE, reconstructionAccuracyE, reconstructionTimeE, reconstructionsE = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"Entropy", "SLSQP")
solC, _,_,_,_ = runSystem(systemID,numTimeSteps,timeStepLength,initState,numStates,numMoments4Closure,k,"CME", "SLSQP")
momsC = distr2moms1D(solC,numMoments4Closure)
In [ ]:
print(solC[:,13])
In [ ]:
moments = np.array([1,1,0,1.1])
checkMomentFeasibility(moments, 11)
In [ ]:
plt.plot(solZ.y[0,:], alpha = np.random.rand(94))
In [ ]:
np.power(2,3)
In [ ]:
2**2.5
In [ ]:
def directDistrReconstruction(numStates, moments):

    numMoments4Closure = moments.shape[0]

    A = np.zeros((numMoments4Closure+1,numMoments4Closure+1))
    b = np.zeros(numMoments4Closure+1)

    # Normalization Criterion
    b[0] = 0
    A[0,0] = numStates
    A[0,1] = numStates*(numStates-1)/2

    for s in range(2,numMoments4Closure+1):
        for i in range(numStates):
            A[0,s] += np.power(i-moments[0],s)

    # Expected Value Criterion
    b[1] = moments[0]-(numStates-1)/2
    A[1,0] = numStates*(numStates-1)/4
    A[1,1] = numStates*(numStates+1)*(2*numStates+1)/12

    for s in range(2,numMoments4Closure+1):
        for i in range(numStates):
            A[1,s] += i*np.power(i-moments[0],s)/2

    # Higher Moment Criteria

    for r in range(2,numMoments4Closure+1):
        b[r] = moments[r-1]
        for i in range(numStates):
            b[r] -= np.power(i-moments[0],r)/numStates
            A[r,0] += np.power(i-moments[0],r)/2
            A[r,1] += np.power(i-moments[0],r)/2*i

            for s in range(2,numMoments4Closure+1):
                A[r,s] += np.power(i-moments[0],r+s)/2
                
                
    lambdaVec = scipy.linalg.solve(A,b)
    distr = np.ones(numStates)/numStates + lambdaVec[0]/2
    
    for i in range(numStates):
        distr[i] += lambdaVec[1]*i/2
        
        for s in range(2,numMoments4Closure+1):
            distr[i] += np.power(i-moments[0],s)/2*lambdaVec[s]
            
    print("This should be the moments: ")
    print(moments)
    print("This are the moments: ")
    matrixDistr = np.zeros((numStates,1))
    matrixDistr[:,0] = distr
    print(distr2moms1D(matrixDistr,numMoments4Closure))
    
    plt.plot(distr)
    plt.show()
    return distr
In [ ]:
moments = np.array([5,10,15,20])
numStates = 11
l = directDistrReconstruction(numStates, moments)

print(l)
In [ ]:
distr2moms1D(CMESol,s)