<s> r"""<STR_LIT>""" <EOL> from __future__ import print_function <EOL> import numpy as np <EOL> def jacobian ( problem , p = None , step = None ) : <EOL> """<STR_LIT>""" <EOL> p_init = problem . getp ( ) <EOL> if p is None : <EOL> p = p_init <EOL> p = np . asarray ( p ) <EOL> bounds = getattr ( problem , '<STR_LIT>' , lambda : None ) ( ) <EOL> def f ( p ) : <EOL> problem . setp ( p ) <EOL> return problem . residuals ( ) <EOL> J = _jacobian_forward ( f , p , bounds , eps = step ) <EOL> problem . setp ( p_init ) <EOL> return J <EOL> def _jacobian_forward ( f , p , bounds , eps = None ) : <EOL> n = len ( p ) <EOL> step = <NUM_LIT> if eps is None else np . sqrt ( eps ) <EOL> fx = f ( p ) <EOL> h = abs ( p ) * step <EOL> h [ h == <NUM_LIT:0> ] = step <EOL> if bounds is not None : <EOL> h [ h + p > bounds [ <NUM_LIT:1> ] ] *= - <NUM_LIT:1.0> <EOL> ee = np . diag ( h ) <EOL> J = [ ] <EOL> for i in range ( n ) : <EOL> J . append ( ( f ( p + ee [ i , : ] ) - fx ) / h [ i ] ) <EOL> return np . vstack ( J ) . T <EOL> def _jacobian_central ( f , p , bounds , eps = None ) : <EOL> n = len ( p ) <EOL> step = <NUM_LIT> if eps is None else np . sqrt ( eps ) <EOL> h = abs ( p ) * step <EOL> h [ h == <NUM_LIT:0> ] = step <EOL> ee = np . diag ( h ) <EOL> J = [ ] <EOL> for i in range ( n ) : <EOL> J . append ( ( f ( p + ee [ i , : ] ) - f ( p - ee [ i , : ] ) ) / ( <NUM_LIT> * h [ i ] ) ) <EOL> return np . vstack ( J ) . T <EOL> def hessian ( problem , p = None , step = None ) : <EOL> """<STR_LIT>""" <EOL> p_init = problem . getp ( ) <EOL> if p is None : <EOL> p = p_init <EOL> p = np . asarray ( p ) <EOL> bounds = getattr ( problem , '<STR_LIT>' , lambda : None ) ( ) <EOL> H = _hessian_forward ( problem . nllf , p , bounds = bounds , eps = step ) <EOL> problem . setp ( p_init ) <EOL> return H <EOL> def _hessian_forward ( f , p , bounds , eps = None ) : <EOL> """<STR_LIT>""" <EOL> n = len ( p ) <EOL> step = <NUM_LIT> if eps is None else np . sqrt ( eps ) <EOL> fx = f ( p ) <EOL> h = abs ( p ) * step <EOL> h [ h == <NUM_LIT:0> ] = step <EOL> if bounds is not None : <EOL> h [ h + p > bounds [ <NUM_LIT:1> ] ] *= - <NUM_LIT:1.0> <EOL> ee = np . diag ( h ) <EOL> g = np . empty ( n , '<STR_LIT:d>' ) <EOL> for i in range ( n ) : <EOL> g [ i ] = f ( p + ee [ i , : ] ) <EOL> H = np . empty ( ( n , n ) , '<STR_LIT:d>' ) <EOL> for i in range ( n ) : <EOL> for j in range ( i , n ) : <EOL> fx_ij = f ( p + ee [ i , : ] + ee [ j , : ] ) <EOL> H [ i , j ] = ( fx_ij - g [ i ] - g [ j ] + fx ) / ( h [ i ] * h [ j ] ) <EOL> H [ j , i ] = H [ i , j ] <EOL> return H <EOL> def _hessian_central ( f , p , bounds , eps = None ) : <EOL> """<STR_LIT>""" <EOL> n = len ( p ) <EOL> step = <NUM_LIT> if eps is None else np . sqrt ( eps ) <EOL> fx = f ( p ) <EOL> h = abs ( p ) * step <EOL> h [ h == <NUM_LIT:0> ] = step <EOL> ee = np . diag ( h ) <EOL> gp = np . empty ( n , '<STR_LIT:d>' ) <EOL> gm = np . empty ( n , '<STR_LIT:d>' ) <EOL> for i in range ( n ) : <EOL> gp [ i ] = f ( p + ee [ i , : ] ) <EOL> gm [ i ] = f ( p - ee [ i , : ] ) <EOL> H = np . empty ( ( n , n ) , '<STR_LIT:d>' ) <EOL> for i in range ( n ) : <EOL> for j in range ( i , n ) : <EOL> fp_ij = f ( p + ee [ i , : ] + ee [ j , : ] ) <EOL> fm_ij = f ( p - ee [ i , : ] - ee [ j , : ] ) <EOL> H [ i , j ] = ( fp_ij - gp [ i ] - gp [ j ] + fm_ij - gm [ i ] - gm [ j ] + <NUM_LIT> * fx ) / ( <NUM_LIT> * h [ i ] * h [ j ] ) <EOL> H [ j , i ] = H [ i , j ] <EOL> return H <EOL> def perturbed_hessian ( H , scale = None ) : <EOL> """<STR_LIT>""" <EOL> from . quasinewton import modelhess <EOL> n = H . shape [ <NUM_LIT:0> ] <EOL> if scale is None : <EOL> scale = np . ones ( n ) <EOL> macheps = np . finfo ( '<STR_LIT:d>' ) . eps <EOL> return modelhess ( n , scale , macheps , H ) <EOL> def chol_stderr ( L ) : <EOL> """<STR_LIT>""" <EOL> return np . sqrt ( <NUM_LIT:1.> / np . diag ( L ) ) <EOL> def chol_cov ( L ) : <EOL> """<STR_LIT>""" <EOL> Linv = np . linalg . inv ( L ) <EOL> return np . dot ( Linv . T . conj ( ) , Linv ) <EOL> def cov ( J , tol = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> u , s , vh = np . linalg . svd ( J , <NUM_LIT:0> ) <EOL> s [ s <= tol ] = tol <EOL> JTJinv = np . dot ( vh . T . conj ( ) / s ** <NUM_LIT:2> , vh ) <EOL> return JTJinv <EOL> def corr ( C ) : <EOL> """<STR_LIT>""" <EOL> Dinv = <NUM_LIT:1.> / stderr ( cov ) <EOL> return np . dot ( Dinv , np . dot ( cov , Dinv ) ) <EOL> def max_correlation ( Rsq ) : <EOL> """<STR_LIT>""" <EOL> return np . max ( np . tril ( Rsq , k = - <NUM_LIT:1> ) ) <EOL> def stderr ( C ) : <EOL> r"""<STR_LIT>""" <EOL> return np . sqrt ( np . diag ( C ) ) <EOL> def demo_hessian ( ) : <EOL> rosen = lambda x : ( <NUM_LIT:1.> - x [ <NUM_LIT:0> ] ) ** <NUM_LIT:2> + <NUM_LIT> * ( x [ <NUM_LIT:1> ] - x [ <NUM_LIT:0> ] ** <NUM_LIT:2> ) ** <NUM_LIT:2> <EOL> p = np . array ( [ <NUM_LIT:1.> , <NUM_LIT:1.> ] ) <EOL> H = _hessian_forward ( rosen , p , bounds = None , eps = <NUM_LIT> ) <EOL> print ( "<STR_LIT>" , H ) <EOL> H = _hessian_central ( rosen , p , bounds = None , eps = <NUM_LIT> ) <EOL> print ( "<STR_LIT>" , H ) <EOL> def demo_jacobian ( ) : <EOL> y = np . array ( [ <NUM_LIT:1.> , <NUM_LIT> , <NUM_LIT> ] ) <EOL> f = lambda x : x [ <NUM_LIT:0> ] * y + x [ <NUM_LIT:1> ] <EOL> p = np . array ( [ <NUM_LIT> , <NUM_LIT> ] ) <EOL> J = _jacobian_forward ( f , p , bounds = None , eps = <NUM_LIT> ) <EOL> print ( "<STR_LIT>" , J ) <EOL> J = _jacobian_central ( f , p , bounds = None , eps = <NUM_LIT> ) <EOL> print ( "<STR_LIT>" , J ) <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> demo_hessian ( ) <EOL> demo_jacobian ( ) </s>
