<s> from __future__ import division <EOL> from __future__ import print_function <EOL> from past . builtins import cmp <EOL> from past . utils import old_div <EOL> import numpy as np <EOL> import matplotlib . pyplot as plt <EOL> import sys <EOL> from math import sqrt <EOL> import scipy . linalg . lapack as lapack <EOL> def jitchol ( A , maxtries = <NUM_LIT:5> ) : <EOL> '''<STR_LIT>''' <EOL> A = np . asfortranarray ( A ) <EOL> L , info = lapack . dpotrf ( A , lower = <NUM_LIT:1> ) <EOL> if info == <NUM_LIT:0> : <EOL> return L <EOL> else : <EOL> diagA = np . diag ( A ) <EOL> if np . any ( diagA <= <NUM_LIT:0.> ) : <EOL> raise np . linalg . LinAlgError ( "<STR_LIT>" ) <EOL> jitter = diagA . mean ( ) * <NUM_LIT> <EOL> while maxtries > <NUM_LIT:0> and np . isfinite ( jitter ) : <EOL> print ( '<STR_LIT>' . format ( jitter ) ) <EOL> try : <EOL> return np . linalg . cholesky ( A + np . eye ( A . shape [ <NUM_LIT:0> ] ) . T * jitter , lower = True ) <EOL> except : <EOL> jitter *= <NUM_LIT:10> <EOL> finally : <EOL> maxtries -= <NUM_LIT:1> <EOL> raise np . linalg . LinAlgError ( "<STR_LIT>" ) <EOL> def solve_chol ( L , B ) : <EOL> '''<STR_LIT>''' <EOL> try : <EOL> assert ( L . shape [ <NUM_LIT:0> ] == L . shape [ <NUM_LIT:1> ] and L . shape [ <NUM_LIT:0> ] == B . shape [ <NUM_LIT:0> ] ) <EOL> except AssertionError : <EOL> raise Exception ( '<STR_LIT>' ) ; <EOL> X = np . linalg . solve ( L , np . linalg . solve ( L . T , B ) ) <EOL> return X <EOL> def unique ( x ) : <EOL> '''<STR_LIT>''' <EOL> y = [ item for sublist in x for item in sublist ] <EOL> if isinstance ( x , np . ndarray ) : <EOL> n , D = x . shape <EOL> assert ( D == <NUM_LIT:1> ) <EOL> y = np . array ( list ( set ( x [ : , <NUM_LIT:0> ] ) ) ) <EOL> y = np . reshape ( y , ( len ( y ) , <NUM_LIT:1> ) ) <EOL> else : <EOL> y = list ( set ( y ) ) <EOL> return y <EOL> def brentmin ( xlow , xupp , Nitmax , tol , f , nout = None , * args ) : <EOL> '''<STR_LIT>''' <EOL> if nout == None : <EOL> nout = <NUM_LIT:0> <EOL> eps = sys . float_info . epsilon <EOL> tol = max ( tol , eps ) <EOL> vargout = f ( xlow , * args ) ; fa = vargout [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] <EOL> vargout = f ( xupp , * args ) ; fb = vargout [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] <EOL> funccount = <NUM_LIT:2> <EOL> seps = sqrt ( eps ) <EOL> c = <NUM_LIT:0.5> * ( <NUM_LIT> - sqrt ( <NUM_LIT> ) ) <EOL> a = xlow <EOL> b = xupp <EOL> v = a + c * ( b - a ) <EOL> w = v <EOL> xf = v <EOL> d = <NUM_LIT:0.> <EOL> e = <NUM_LIT:0.> <EOL> x = xf <EOL> vargout = f ( x , * args ) <EOL> fx = vargout [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] <EOL> varargout = vargout [ <NUM_LIT:1> : ] <EOL> funccount += <NUM_LIT:1> <EOL> fv = fx ; fw = fx <EOL> xm = <NUM_LIT:0.5> * ( a + b ) <EOL> tol1 = seps * abs ( xf ) + old_div ( tol , <NUM_LIT> ) ; <EOL> tol2 = <NUM_LIT> * tol1 <EOL> while ( abs ( xf - xm ) > ( tol2 - <NUM_LIT:0.5> * ( b - a ) ) ) : <EOL> gs = True <EOL> if abs ( e ) > tol1 : <EOL> gs = False <EOL> r = ( xf - w ) * ( fx - fv ) <EOL> q = ( xf - v ) * ( fx - fw ) <EOL> p = ( xf - v ) * q - ( xf - w ) * r <EOL> q = <NUM_LIT> * ( q - r ) <EOL> if q > <NUM_LIT:0.0> : <EOL> p = - p <EOL> q = abs ( q ) <EOL> r = e ; e = d <EOL> if ( ( abs ( p ) < abs ( <NUM_LIT:0.5> * q * r ) ) and ( p > q * ( a - xf ) ) and ( p < q * ( b - xf ) ) ) : <EOL> d = old_div ( p , q ) <EOL> x = xf + d <EOL> if ( ( x - a ) < tol2 ) or ( ( b - x ) < tol2 ) : <EOL> si = cmp ( xm - xf , <NUM_LIT:0> ) <EOL> if ( ( xm - xf ) == <NUM_LIT:0> ) : si += <NUM_LIT:1> <EOL> d = tol1 * si <EOL> else : <EOL> gs = True <EOL> if gs : <EOL> if xf >= xm : e = a - xf <EOL> else : <EOL> e = b - xf <EOL> d = c * e <EOL> si = cmp ( d , <NUM_LIT:0> ) <EOL> if ( d == <NUM_LIT:0> ) : si += <NUM_LIT:1> <EOL> x = xf + si * max ( abs ( d ) , tol1 ) <EOL> vargout = f ( x , * args ) ; fu = vargout [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] ; varargout = vargout [ <NUM_LIT:1> : ] <EOL> funccount += <NUM_LIT:1> <EOL> if fu <= fx : <EOL> if x >= xf : a = xf <EOL> else : b = xf <EOL> v = w ; fv = fw <EOL> w = xf ; fw = fx <EOL> xf = x ; fx = fu <EOL> else : <EOL> if x < xf : <EOL> a = x <EOL> else : <EOL> b = x <EOL> if ( ( fu <= fw ) or ( w == xf ) ) : <EOL> v = w ; fv = fw <EOL> w = x ; fw = fu <EOL> elif ( ( fu <= fv ) or ( ( v == xf ) or ( v == w ) ) ) : <EOL> v = x ; fv = fu <EOL> xm = <NUM_LIT:0.5> * ( a + b ) <EOL> tol1 = seps * abs ( xf ) + old_div ( tol , <NUM_LIT> ) ; tol2 = <NUM_LIT> * tol1 <EOL> if funccount >= Nitmax : <EOL> break <EOL> if ( ( fa < fx ) and ( fa <= fb ) ) : <EOL> xf = xlow ; fx = fa <EOL> elif fb < fx : <EOL> xf = xupp ; fx = fb <EOL> fmin = fx <EOL> xmin = xf <EOL> vargout = [ xmin , fmin , funccount ] <EOL> for vv in varargout : <EOL> vargout . append ( vv ) <EOL> return vargout <EOL> def cholupdate ( R , x , sgn = '<STR_LIT:+>' ) : <EOL> '''<STR_LIT>''' <EOL> if len ( x . shape ) == <NUM_LIT:1> : <EOL> x = np . reshape ( x , ( x . shape [ <NUM_LIT:0> ] , <NUM_LIT:1> ) ) <EOL> assert ( R . shape [ <NUM_LIT:0> ] == x . shape [ <NUM_LIT:0> ] ) <EOL> A = np . dot ( R . T , R ) <EOL> if sgn == '<STR_LIT:+>' : <EOL> R1 = A + np . dot ( x , x . T ) <EOL> elif sgn == '<STR_LIT:->' : <EOL> R1 = A - np . dot ( x , x . T ) <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> return jitchol ( R1 ) . T </s>
