python_header.py

import numpy
from numpy import linalg as LA
from numpy.random import random
from scipy.optimize import curve_fit
from scipy.spatial import cKDTree
import matplotlib as mpl
import matplotlib.patches as patches
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import seaborn as sns
import pylab as pylab
import glob 
from glob import glob
import matplotlib.ticker as mtick
import matplotlib.colors as mcol
import time, math
import sys, os
from scipy import optimize
from matplotlib.ticker import AutoMinorLocator
import matplotlib.cm as cm
import pandas as pd

def SolveC(A, M, mass_inv, b, maxiter):
    d = mass_inv @ b
    x = d
    normr = []

    for i in range(1, maxiter):
        d = A @ d
        x += d

        r = M @ x - b
        normr.append((numpy.dot(r,r)))

        #print("iteration {}, ||r||^2={:.3e}".format(i, normr[-1]))
    return x,normr


def LinearCG(A, b, x0, tol=1e-5, maxiter=10):
    xk = x0
    rk = numpy.dot(A, xk) - b
    pk = -rk
    rk_norm = numpy.linalg.norm(rk)
    residuals = [rk_norm]
    
    num_iter = 0
    curve_x = [xk]
    while rk_norm > tol:
        apk = numpy.dot(A, pk)
        rkrk = numpy.dot(rk, rk)
        
        alpha = rkrk / numpy.dot(pk, apk)
        xk = xk + alpha * pk
        rk = rk + alpha * apk
        beta = numpy.dot(rk, rk) / rkrk
        pk = -rk + beta * pk
        
        num_iter += 1
        curve_x.append(xk)
        rk_norm = numpy.linalg.norm(rk)
        residuals.append(rk_norm)
        #print('Iteration: {} \t residual = {:.2e}'.
        #      format(num_iter, rk_norm))
        if num_iter >= maxiter:
            break
    
    #print('\nSolution: \t x = {}'.format(xk))
        
    return xk,residuals

def PreconditionedLinearCG(A, b, invM, tol=1e-5, maxiter=10):
    xk = invM @ b
    rk = b - numpy.dot(A, xk)
    pk = invM @ rk
    rk_norm = numpy.linalg.norm(rk)
    residuals = [numpy.dot(rk, rk)]
    
    num_iter = 0
    curve_x = [xk]
    while rk_norm > tol:
        apk = numpy.dot(A, pk)
        rkrk = numpy.dot(rk, invM @ rk)
        
        alpha = rkrk / numpy.dot(pk, apk)
        xk = xk + alpha * pk
        rk = rk - alpha * apk
        beta = numpy.dot(rk, invM @ rk) / rkrk
        pk = invM @ rk + beta * pk
        
        num_iter += 1
        curve_x.append(xk)
        rk_norm = numpy.linalg.norm(rk)
        R = b - A @ xk
        residuals.append(numpy.dot(R,R))
        #print('Iteration: {} \t residual = {:.2e}'.
        #      format(num_iter, rk_norm))
        if num_iter >= maxiter:
            break
    
    #print('\nSolution: \t x = {}'.format(xk))
        
    return xk,residuals

def PreconditionedLinearCG2(A, b, invM, tol=1e-5, maxiter=10):
    xk = b
    rk = invM @ (b - A @ (invM @ xk))
    pk = rk
    rk_norm = numpy.linalg.norm(rk)
    R = b - A @ (invM @ xk)
    residuals = [numpy.dot(R,R)]
    
    num_iter = 0
    curve_x = [xk]
    #print('Iteration: {} \t residual = {:.2e}'.
    #      format(num_iter, residuals[-1]))
    while rk_norm > tol:
        apk = invM @ (A @ (invM @ pk))
        rkrk = numpy.dot(rk, rk)

        pk_dot_apk = numpy.dot(pk, apk)
        alpha = rkrk / numpy.dot(pk, apk)
        xk = xk + alpha * pk
        rk = rk - alpha * apk
        beta = numpy.dot(rk, rk) / rkrk
        pk = rk + beta * pk
        
        num_iter += 1
        curve_x.append(xk)
        rk_norm = numpy.linalg.norm(rk)
        R = b - A @ (invM @ xk)
        residuals.append(numpy.dot(R,R))
        #print('Iteration: {} \t alpha = {} \t rkrk={:.5e} \t pk_dot_apk={:.5e} \t residual = {:.2e}'.
        #      format(num_iter, alpha, rkrk, pk_dot_apk, residuals[-1]))
        if num_iter >= maxiter:
            break
    
    #print('\nSolution: \t x = {}'.format(xk))
        
    return invM @ xk,residuals

def LinearCG2(A, b, x0, tol=1e-5, maxiter=10):
    xk = x0
    rk = numpy.dot(A, xk) - b
    pk = -rk
    rk_norm = numpy.linalg.norm(rk)
    
    num_iter = 0
    curve_x = [xk]
    residuals = []
    i = 0
    while rk_norm > tol:
        apk = numpy.dot(A, pk)
        rkrk = numpy.dot(rk, rk)
        
        alpha = rkrk / numpy.dot(pk, apk)
        xk = xk + alpha * pk
        if i < 10:
            rk = rk + alpha * apk
            i += 1
        else:
            rk = numpy.dot(A, xk) - b
            i = 0
        beta = numpy.dot(rk, rk) / rkrk
        pk = -rk + beta * pk
        
        num_iter += 1
        curve_x.append(xk)
        rk_norm = numpy.linalg.norm(rk)
        residuals.append(rk_norm)
        #print('Iteration: {} \t residual = {:.2e}'.
        #      format(num_iter, rk_norm))
        if num_iter >= maxiter:
            break
    
    #print('\nSolution: \t x = {}'.format(xk))
        
    return xk,residuals

def find_cells_per_unit_length(f):
    posI = f.find('N_')
    posE = f.find('/log.mpm')
    return int(f[posI+2:posE])

def update_gitignore(pdfname):
    pdfline = "/" + pdfname 
    present = False 
    with open('.gitignore', 'r+') as f: 
        for line in f: 
            if pdfline in line: 
                present = True 
                break 
        if  not present: 
            f.write(pdfline + "\n") 

def label_line(ax, line, label, color='0.5', fs=14, halign='left'):
    """Add an annotation to the given line with appropriate placement and rotation.
    Based on code from:
        [How to rotate matplotlib annotation to match a line?]
        (http://stackoverflow.com/a/18800233/230468)
        User: [Adam](http://stackoverflow.com/users/321772/adam)
    Arguments
    ---------
    ax : `matplotlib.axes.Axes` object
        Axes on which the label should be added.
    line : `matplotlib.lines.Line2D` object
        Line which is being labelled.
    label : str
        Text which should be drawn as the label.
    ...
    Returns
    -------
    text : `matplotlib.text.Text` object
    """
    xdata, ydata = line.get_data()
    x1 = xdata[0]
    x2 = xdata[-1]
    y1 = ydata[0]
    y2 = ydata[-1]

    if halign.startswith('l'):
        xx = x1
        halign = 'left'
    elif halign.startswith('r'):
        xx = x2
        halign = 'right'
    elif halign.startswith('c'):
        xx = 0.5*(x1 + x2)
        halign = 'center'
    else:
        raise ValueError("Unrecogznied `halign` = '{}'.".format(halign))

    yy = numpy.interp(xx, xdata, ydata)

    ylim = ax.get_ylim()
    # xytext = (10, 10)
    xytext = (0, 0)
    text = ax.annotate(label, xy=(xx, yy), xytext=xytext, textcoords='offset points',
                       size=fs, color=color, zorder=1,
                       horizontalalignment=halign, verticalalignment='center_baseline')

    sp1 = ax.transData.transform_point((x1, y1))
    sp2 = ax.transData.transform_point((x2, y2))

    rise = (sp2[1] - sp1[1])
    run = (sp2[0] - sp1[0])

    slope_degrees = numpy.degrees(numpy.arctan2(rise, run))
    text.set_rotation_mode('anchor')
    text.set_rotation(slope_degrees)
    ax.set_ylim(ylim)
    return text

def set_size(width, fraction=1):
    """ Set aesthetic figure dimensions to avoid scaling in latex.

    Parameters
    ----------
    width: float
            Width in pts
    fraction: float
            Fraction of the width which you wish the figure to occupy

    Returns
    -------
    fig_dim: tuple
            Dimensions of figure in inches
    """

    if width == 'elsevier':
        width_pt = 468.
    elif width == 'beamer':
        width_pt = 307.28987
    else:
        width_pt = width

    # Width of figure
    fig_width_pt = width_pt * fraction

    # Convert from pt to inches
    inches_per_pt = 1 / 72.27

    # Golden ratio to set aesthetic figure height
    golden_ratio = 0.75 # (5**.5 - 1) / 2

    # Figure width in inches
    fig_width_in = fig_width_pt * inches_per_pt
    # Figure height in inches
    fig_height_in = fig_width_in * golden_ratio

    fig_dim = (fig_width_in, fig_height_in)

    return fig_dim

def power_law(x, C, alpha):
    return C * numpy.power(x, -alpha)

# Using seaborn's style
#plt.style.use('seaborn-colorblind') # dark_background, 
sns.set_style('ticks') # dark_background, 

width = 345

nice_fonts = {
    # Use LaTeX to write all text
    "text.usetex": True,
    #"font.family": "avant",
    # Use 10pt font in plots, to match 10pt font in document
    "axes.labelsize": 8,
    "font.size": 8,
    # Make the legend/label fonts a little smaller
    "legend.fontsize": 8,
    "xtick.labelsize": 8,
    "ytick.labelsize": 8,
}

plt.rcParams.update(nice_fonts)

markersize = 3