signature:
def circles ( st , layer , axis , ax = None , talpha = 1.0 , cedge = 'white' , cface = 'white' ) :
    """Plots a set of circles corresponding to a slice through the platonic
 structure . Copied from twoslice _ overlay with comments , standaloneness .
 Inputs
 pos : array of particle positions ; [ N , 3]
 rad : array of particle radii ; [ N ]
 ax : plt . axis instance
 layer : Which layer of the slice to use .
 axis : The slice of the image , 0 , 1 , or 2.
 cedge : edge color
 cface : face color
 talpha : Alpha of the thing"""

implementation:
pos = st . obj_get_positions ( )
rad = st . obj_get_radii ( )
shape = st . ishape . shape . tolist ( )
shape . pop ( axis )
# shape is now the shape of the image
if ax is None :
    fig = plt . figure ( )
    axisbg = 'white' if cface == 'black' else 'black'
    sx , sy = ( ( 1 , shape [ 1 ] / float ( shape [ 0 ] ) ) if shape [ 0 ] > shape [ 1 ] else ( shape [ 0 ] / float ( shape [ 1 ] ) , 1 ) )
    ax = fig . add_axes ( ( 0 , 0 , sx , sy ) , axisbg = axisbg )
# get the index of the particles we want to include
particles = np . arange ( len ( pos ) ) [ np . abs ( pos [ : , axis ] - layer ) < rad ]
# for each of these particles display the effective radius
# in the proper place
scale = 1.0
# np . max ( shape ) . astype ( ' float ' )
for i in particles :
    p = pos [ i ] . copy ( )
    r = 2 * np . sqrt ( rad [ i ] ** 2 - ( p [ axis ] - layer ) ** 2 )
    # CIRCLE IS IN FIGURE COORDINATES ! ! !
    if axis == 0 :
        ix = 1 ;
        iy = 2
    elif axis == 1 :
        ix = 0 ;
        iy = 2
    elif axis == 2 :
        ix = 0 ;
        iy = 1
    c = Circle ( ( p [ ix ] / scale , p [ iy ] / scale ) , radius = r / 2 / scale , fc = cface , ec = cedge , alpha = talpha )
    ax . add_patch ( c )
# plt . axis ( [ 0,1,0,1 ] )
plt . axis ( 'equal' )
# circles not ellipses
return ax


signature:
def missing_particle ( separation = 0.0 , radius = RADIUS , SNR = 20 ) :
    """create a two particle state and compare it to featuring using a single particle guess"""

implementation:
# create a base image of one particle
s = init . create_two_particle_state ( imsize = 6 * radius + 4 , axis = 'x' , sigma = 1.0 / SNR , delta = separation , radius = radius , stateargs = { 'varyn' : True } , psfargs = { 'error' : 1e-6 } )
s . obj . typ [ 1 ] = 0.
s . reset ( )
return s , s . obj . pos . copy ( )


signature:
def get_rand_Japprox ( s , params , num_inds = 1000 , include_cost = False , ** kwargs ) :
    """Calculates a random approximation to J by returning J only at a
 set of random pixel / voxel locations .
 Parameters
 s : : class : ` peri . states . State `
 The state to calculate J for .
 params : List
 The list of parameter names to calculate the gradient of .
 num _ inds : Int , optional .
 The number of pix / voxels at which to calculate the random
 approximation to J . Default is 1000.
 include _ cost : Bool , optional
 Set to True to append a finite - difference measure of the full
 cost gradient onto the returned J .
 Other Parameters
 All kwargs parameters get passed to s . gradmodel only .
 Returns
 J : numpy . ndarray
 [ d , num _ inds ] array of J , at the given indices .
 return _ inds : numpy . ndarray or slice
 [ num _ inds ] element array or slice ( 0 , None ) of the model
 indices at which J was evaluated ."""

implementation:
start_time = time . time ( )
tot_pix = s . residuals . size
if num_inds < tot_pix :
    inds = np . random . choice ( tot_pix , size = num_inds , replace = False )
    slicer = None
    return_inds = np . sort ( inds )
else :
    inds = None
    return_inds = slice ( 0 , None )
    slicer = [ slice ( 0 , None ) ] * len ( s . residuals . shape )
if include_cost :
    Jact , ge = s . gradmodel_e ( params = params , inds = inds , slicer = slicer , flat = False , ** kwargs )
    Jact *= - 1
    J = [ Jact , ge ]
else :
    J = - s . gradmodel ( params = params , inds = inds , slicer = slicer , flat = False , ** kwargs )
CLOG . debug ( 'J:\t%f' % ( time . time ( ) - start_time ) )
return J , return_inds


signature:
def name_globals ( s , remove_params = None ) :
    """Returns a list of the global parameter names .
 Parameters
 s : : class : ` peri . states . ImageState `
 The state to name the globals of .
 remove _ params : Set or None
 A set of unique additional parameters to remove from the globals
 list .
 Returns
 all _ params : list
 The list of the global parameter names , with each of
 remove _ params removed ."""

implementation:
all_params = s . params
for p in s . param_particle ( np . arange ( s . obj_get_positions ( ) . shape [ 0 ] ) ) :
    all_params . remove ( p )
if remove_params is not None :
    for p in set ( remove_params ) :
        all_params . remove ( p )
return all_params


signature:
def get_num_px_jtj ( s , nparams , decimate = 1 , max_mem = 1e9 , min_redundant = 20 ) :
    """Calculates the number of pixels to use for J at a given memory usage .
 Tries to pick a number of pixels as ( size of image / ` decimate ` ) .
 However , clips this to a maximum size and minimum size to ensure that
 (1 ) too much memory isn ' t used and ( 2 ) J has enough elements so that
 the inverse of JTJ will be well - conditioned .
 Parameters
 s : : class : ` peri . states . State `
 The state on which to calculate J .
 nparams : Int
 The number of parameters that will be included in J .
 decimate : Int , optional
 The amount to decimate the number of pixels in the image by ,
 i . e . tries to pick num _ px = size of image / decimate .
 Default is 1
 max _ mem : Numeric , optional
 The maximum allowed memory , in bytes , for J to occupy at
 double - precision . Default is 1e9.
 min _ redundant : Int , optional
 The number of pixels must be at least ` min _ redundant ` *
 ` nparams ` . If not , an error is raised . Default is 20
 Returns
 num _ px : Int
 The number of pixels at which to calcualte J ."""

implementation:
# 1 . Max for a given max _ mem :
px_mem = int ( max_mem // 8 // nparams )
# 1 float = 8 bytes
# 2 . num _ pix for a given redundancy
px_red = min_redundant * nparams
# 3 . And # desired for decimation
px_dec = s . residuals . size // decimate
if px_red > px_mem :
    raise RuntimeError ( 'Insufficient max_mem for desired redundancy.' )
num_px = np . clip ( px_dec , px_red , px_mem ) . astype ( 'int' )
return num_px


signature:
def vectorize_damping ( params , damping = 1.0 , increase_list = [ [ 'psf-' , 1e4 ] ] ) :
    """Returns a non - constant damping vector , allowing certain parameters to be
 more strongly damped than others .
 Parameters
 params : List
 The list of parameter names , in order .
 damping : Float
 The default value of the damping .
 increase _ list : List
 A nested 2 - element list of the params to increase and their
 scale factors . All parameters containing the string
 increase _ list [ i ] [ 0 ] are increased by a factor increase _ list [ i ] [ 1 ] .
 Returns
 damp _ vec : np . ndarray
 The damping vector to use ."""

implementation:
damp_vec = np . ones ( len ( params ) ) * damping
for nm , fctr in increase_list :
    for a in range ( damp_vec . size ) :
        if nm in params [ a ] :
            damp_vec [ a ] *= fctr
return damp_vec


signature:
def low_mem_sq ( m , step = 100000 ) :
    """np . dot ( m , m . T ) with low mem usage , by doing it in small steps"""

implementation:
if not m . flags . c_contiguous :
    raise ValueError ( 'm must be C ordered for this to work with less mem.' )
# - - can make this even faster with pre - allocating arrays , but not worth it
# right now
# mmt = np . zeros ( [ m . shape [ 0 ] , m . shape [ 0 ] ] ) # 6us
# mt _ tmp = np . zeros ( [ step , m . shape [ 0 ] ] )
# for a in range ( 0 , m . shape [ 1 ] , step ) :
# mx = min ( a + step , m . shape [ 1 ] )
# mt _ tmp [ : mx - a , : ] = m . T [ a : mx ]
# # np . dot ( m _ tmp , m . T , out = mmt [ a : mx ] )
# # np . dot ( m , m [ a : mx ] . T , out = mmt [ : , a : mx ] )
# np . dot ( m [ : , a : mx ] , mt _ tmp [ : mx ] , out = mmt )
# return mmt
mmt = np . zeros ( [ m . shape [ 0 ] , m . shape [ 0 ] ] )
# 6us
# m _ tmp = np . zeros ( [ step , m . shape [ 1 ] ] )
for a in range ( 0 , m . shape [ 0 ] , step ) :
    mx = min ( a + step , m . shape [ 1 ] )
    # m _ tmp [ : ] = m [ a : mx ]
    # np . dot ( m _ tmp , m . T , out = mmt [ a : mx ] )
    mmt [ : , a : mx ] = np . dot ( m , m [ a : mx ] . T )
return mmt


signature:
def find_particles_in_tile ( positions , tile ) :
    """Finds the particles in a tile , as numpy . ndarray of ints .
 Parameters
 positions : ` numpy . ndarray `
 [ N , 3 ] array of the particle positions to check in the tile
 tile : : class : ` peri . util . Tile ` instance
 Tile of the region inside which to check for particles .
 Returns
 numpy . ndarray , int
 The indices of the particles in the tile ."""

implementation:
bools = tile . contains ( positions )
return np . arange ( bools . size ) [ bools ]


signature:
def separate_particles_into_groups ( s , region_size = 40 , bounds = None , doshift = False ) :
    """Separates particles into convenient groups for optimization .
 Given a state , returns a list of groups of particles . Each group of
 particles are located near each other in the image . Every particle
 located in the desired region is contained in exactly 1 group .
 Parameters
 s : : class : ` peri . states . ImageState `
 The peri state to find particles in .
 region _ size : Int or 3 - element list - like of ints , optional
 The size of the box . Groups particles into boxes of shape
 ( region _ size [ 0 ] , region _ size [ 1 ] , region _ size [ 2 ] ) . If region _ size
 is a scalar , the box is a cube of length region _ size .
 Default is 40.
 bounds : 2 - element list - like of 3 - element lists , optional
 The sub - region of the image over which to look for particles .
 bounds [ 0 ] : The lower - left corner of the image region .
 bounds [ 1 ] : The upper - right corner of the image region .
 Default ( None - > ( [ 0,0,0 ] , s . oshape . shape ) ) is a box of the entire
 image size , i . e . the default places every particle in the image
 somewhere in the groups .
 doshift : { True , False , ` ' rand ' ` } , optional
 Whether or not to shift the tile boxes by half a region size , to
 prevent the same particles to be chosen every time . If ` ' rand ' ` ,
 randomly chooses either True or False . Default is False
 Returns
 particle _ groups : List
 Each element of particle _ groups is an int numpy . ndarray of the
 group of nearby particles . Only contains groups with a nonzero
 number of particles , so the elements don ' t necessarily correspond
 to a given image region ."""

implementation:
imtile = s . oshape . translate ( - s . pad )
bounding_tile = ( imtile if bounds is None else Tile ( bounds [ 0 ] , bounds [ 1 ] ) )
rs = ( np . ones ( bounding_tile . dim , dtype = 'int' ) * region_size if np . size ( region_size ) == 1 else np . array ( region_size ) )
n_translate = np . ceil ( bounding_tile . shape . astype ( 'float' ) / rs ) . astype ( 'int' )
particle_groups = [ ]
tile = Tile ( left = bounding_tile . l , right = bounding_tile . l + rs )
if doshift == 'rand' :
    doshift = np . random . choice ( [ True , False ] )
if doshift :
    shift = rs // 2
    n_translate += 1
else :
    shift = 0
deltas = np . meshgrid ( * [ np . arange ( i ) for i in n_translate ] )
positions = s . obj_get_positions ( )
if bounds is None : # FIXME this ( deliberately ) masks a problem where optimization
# places particles outside the image . However , it ensures that
# all particles are in at least one group when ` bounds is None ` ,
# which is the use case within opt . The 1e - 3 is to ensure that
# they are inside the box and not on the edge .
    positions = np . clip ( positions , imtile . l + 1e-3 , imtile . r - 1e-3 )
groups = list ( map ( lambda * args : find_particles_in_tile ( positions , tile . translate ( np . array ( args ) * rs - shift ) ) , * [ d . ravel ( ) for d in deltas ] ) )
for i in range ( len ( groups ) - 1 , - 1 , - 1 ) :
    if groups [ i ] . size == 0 :
        groups . pop ( i )
assert _check_groups ( s , groups )
return groups


signature:
def _check_groups ( s , groups ) :
    """Ensures that all particles are included in exactly 1 group"""

implementation:
ans = [ ]
for g in groups :
    ans . extend ( g )
if np . unique ( ans ) . size != np . size ( ans ) :
    return False
elif np . unique ( ans ) . size != s . obj_get_positions ( ) . shape [ 0 ] :
    return False
else :
    return ( np . arange ( s . obj_get_radii ( ) . size ) == np . sort ( ans ) ) . all ( )


signature:
def calc_particle_group_region_size ( s , region_size = 40 , max_mem = 1e9 , ** kwargs ) :
    """Finds the biggest region size for LM particle optimization with a
 given memory constraint .
 Input Parameters
 s : : class : ` peri . states . ImageState `
 The state with the particles
 region _ size : Int or 3 - element list - like of ints , optional .
 The initial guess for the region size . Default is 40
 max _ mem : Numeric , optional
 The maximum memory for the optimizer to take . Default is 1e9
 Other Parameters
 bounds : 2 - element list - like of 3 - element lists .
 The sub - region of the image over which to look for particles .
 bounds [ 0 ] : The lower - left corner of the image region .
 bounds [ 1 ] : The upper - right corner of the image region .
 Default ( None - > ( [ 0,0,0 ] , s . oshape . shape ) ) is a box of the entire
 image size , i . e . the default places every particle in the image
 somewhere in the groups .
 Returns
 region _ size : numpy . ndarray of ints of the region size ."""

implementation:
region_size = np . array ( region_size ) . astype ( 'int' )
def calc_mem_usage ( region_size ) :
    rs = np . array ( region_size )
    particle_groups = separate_particles_into_groups ( s , region_size = rs . tolist ( ) , ** kwargs )
    # The actual mem usage is the max of the memory usage of all the
    # particle groups . However this is too slow . So instead we use the
    # max of the memory of the biggest 5 particle groups :
    numpart = [ np . size ( g ) for g in particle_groups ]
    biggroups = [ particle_groups [ i ] for i in np . argsort ( numpart ) [ - 5 : ] ]
    def get_tile_jsize ( group ) :
        nms = s . param_particle ( group )
        tile = s . get_update_io_tiles ( nms , s . get_values ( nms ) ) [ 2 ]
        return tile . shape . prod ( ) * len ( nms )
    mems = [ 8 * get_tile_jsize ( g ) for g in biggroups ]
    # 8 for bytes / float64
    return np . max ( mems )
im_shape = s . oshape . shape
if calc_mem_usage ( region_size ) > max_mem :
    while ( ( calc_mem_usage ( region_size ) > max_mem ) and np . any ( region_size > 2 ) ) :
        region_size = np . clip ( region_size - 1 , 2 , im_shape )
else :
    while ( ( calc_mem_usage ( region_size ) < max_mem ) and np . any ( region_size < im_shape ) ) :
        region_size = np . clip ( region_size + 1 , 2 , im_shape )
    region_size -= 1
# need to be < memory , so we undo 1 iteration
return region_size


signature:
def get_residuals_update_tile ( st , padded_tile ) :
    """Translates a tile in the padded image to the unpadded image .
 Given a state and a tile that corresponds to the padded image , returns
 a tile that corresponds to the the corresponding pixels of the difference
 image
 Parameters
 st : : class : ` peri . states . State `
 The state
 padded _ tile : : class : ` peri . util . Tile `
 The tile in the padded image .
 Returns
 : class : ` peri . util . Tile `
 The tile corresponding to padded _ tile in the unpadded image ."""

implementation:
inner_tile = st . ishape . intersection ( [ st . ishape , padded_tile ] )
return inner_tile . translate ( - st . pad )


signature:
def find_best_step ( err_vals ) :
    """Returns the index of the lowest of the passed values . Catches nans etc ."""

implementation:
if np . all ( np . isnan ( err_vals ) ) :
    raise ValueError ( 'All err_vals are nans!' )
return np . nanargmin ( err_vals )


signature:
def do_levmarq ( s , param_names , damping = 0.1 , decrease_damp_factor = 10. , run_length = 6 , eig_update = True , collect_stats = False , rz_order = 0 , run_type = 2 , ** kwargs ) :
    """Runs Levenberg - Marquardt optimization on a state .
 Convenience wrapper for LMGlobals . Same keyword args , but the defaults
 have been set to useful values for optimizing globals .
 See LMGlobals and LMEngine for documentation .
 See Also
 do _ levmarq _ particles : Levenberg - Marquardt optimization of a
 specified set of particles .
 do _ levmarq _ all _ particle _ groups : Levenberg - Marquardt optimization
 of all the particles in the state .
 LMGlobals : Optimizer object ; the workhorse of do _ levmarq .
 LMEngine : Engine superclass for all the optimizers ."""

implementation:
if rz_order > 0 :
    aug = AugmentedState ( s , param_names , rz_order = rz_order )
    lm = LMAugmentedState ( aug , damping = damping , run_length = run_length , decrease_damp_factor = decrease_damp_factor , eig_update = eig_update , ** kwargs )
else :
    lm = LMGlobals ( s , param_names , damping = damping , run_length = run_length , decrease_damp_factor = decrease_damp_factor , eig_update = eig_update , ** kwargs )
if run_type == 2 :
    lm . do_run_2 ( )
elif run_type == 1 :
    lm . do_run_1 ( )
else :
    raise ValueError ( 'run_type=1,2 only' )
if collect_stats :
    return lm . get_termination_stats ( )


signature:
def do_levmarq_particles ( s , particles , damping = 1.0 , decrease_damp_factor = 10. , run_length = 4 , collect_stats = False , max_iter = 2 , ** kwargs ) :
    """Levenberg - Marquardt optimization on a set of particles .
 Convenience wrapper for LMParticles . Same keyword args , but the
 defaults have been set to useful values for optimizing particles .
 See LMParticles and LMEngine for documentation .
 See Also
 do _ levmarq _ all _ particle _ groups : Levenberg - Marquardt optimization
 of all the particles in the state .
 do _ levmarq : Levenberg - Marquardt optimization of the entire state ;
 useful for optimizing global parameters .
 LMParticles : Optimizer object ; the workhorse of do _ levmarq _ particles .
 LMEngine : Engine superclass for all the optimizers ."""

implementation:
lp = LMParticles ( s , particles , damping = damping , run_length = run_length , decrease_damp_factor = decrease_damp_factor , max_iter = max_iter , ** kwargs )
lp . do_run_2 ( )
if collect_stats :
    return lp . get_termination_stats ( )


signature:
def do_levmarq_all_particle_groups ( s , region_size = 40 , max_iter = 2 , damping = 1.0 , decrease_damp_factor = 10. , run_length = 4 , collect_stats = False , ** kwargs ) :
    """Levenberg - Marquardt optimization for every particle in the state .
 Convenience wrapper for LMParticleGroupCollection . Same keyword args ,
 but I ' ve set the defaults to what I ' ve found to be useful values for
 optimizing particles . See LMParticleGroupCollection for documentation .
 See Also
 do _ levmarq _ particles : Levenberg - Marquardt optimization of a
 specified set of particles .
 do _ levmarq : Levenberg - Marquardt optimization of the entire state ;
 useful for optimizing global parameters .
 LMParticleGroupCollection : The workhorse of do _ levmarq .
 LMEngine : Engine superclass for all the optimizers ."""

implementation:
lp = LMParticleGroupCollection ( s , region_size = region_size , damping = damping , run_length = run_length , decrease_damp_factor = decrease_damp_factor , get_cos = collect_stats , max_iter = max_iter , ** kwargs )
lp . do_run_2 ( )
if collect_stats :
    return lp . stats


signature:
def do_levmarq_n_directions ( s , directions , max_iter = 2 , run_length = 2 , damping = 1e-3 , collect_stats = False , marquardt_damping = True , ** kwargs ) :
    """Optimization of a state along a specific set of directions in parameter
 space .
 Parameters
 s : : class : ` peri . states . State `
 The state to optimize
 directions : np . ndarray
 [ n , d ] element numpy . ndarray of the n directions in the d -
 dimensional space to optimize along . ` directions ` is trans -
 formed to a unit vector internally
 Other Parameters
 Any parameters passed to LMEngine ."""

implementation:
# normal = direction / np . sqrt ( np . dot ( direction , direction ) )
normals = np . array ( [ d / np . sqrt ( np . dot ( d , d ) ) for d in directions ] )
if np . isnan ( normals ) . any ( ) :
    raise ValueError ( '`directions` must not be 0s or contain nan' )
obj = OptState ( s , normals )
lo = LMOptObj ( obj , max_iter = max_iter , run_length = run_length , damping = damping , marquardt_damping = marquardt_damping , ** kwargs )
lo . do_run_1 ( )
if collect_stats :
    return lo . get_termination_stats ( )


signature:
def burn ( s , n_loop = 6 , collect_stats = False , desc = '' , rz_order = 0 , fractol = 1e-4 , errtol = 1e-2 , mode = 'burn' , max_mem = 1e9 , include_rad = True , do_line_min = 'default' , partial_log = False , dowarn = True ) :
    """Optimizes all the parameters of a state .
 Burns a state through calling LMParticleGroupCollection and LMGlobals /
 LMAugmentedState .
 Parameters
 s : : class : ` peri . states . ImageState `
 The state to optimize
 n _ loop : Int , optional
 The number of times to loop over in the optimizer . Default is 6.
 collect _ stats : Bool , optional
 Whether or not to collect information on the optimizer ' s
 performance . Default is False .
 desc : string , optional
 Description to append to the states . save ( ) call every loop .
 Set to None to avoid saving . Default is ' ' , which selects
 one of ' burning ' , ' polishing ' , ' doing _ positions '
 rz _ order : Int , optional
 Set to an int > 0 to optimize with an augmented state ( R ( z ) as
 a global parameter ) vs . with the normal global parameters ;
 rz _ order is the order of the polynomial approximate for R ( z ) .
 Default is 0 ( no augmented state ) .
 fractol : Float , optional
 Fractional change in error at which to terminate . Default 1e - 4
 errtol : Float , optional
 Absolute change in error at which to terminate . Default 1e - 2
 mode : { ' burn ' , ' do - particles ' , or ' polish ' } , optional
 What mode to optimize with .
 * ' burn ' : Your state is far from the minimum .
 * ' do - particles ' : Positions far from minimum , globals well - fit .
 * ' polish ' : The state is close to the minimum .
 ' burn ' is the default . Only ` polish ` will get to the global
 minimum .
 max _ mem : Numeric , optional
 The maximum amount of memory allowed for the optimizers ' J ' s ,
 for both particles & globals . Default is 1e9 , i . e . 1GB per
 optimizer .
 do _ line _ min : Bool or ' default ' , optional
 Set to True to do an additional , third optimization per loop
 which optimizes along the subspace spanned by the last 3 steps
 of the burn ( ) ' s trajectory . In principle this should signifi -
 cantly speed up the convergence ; in practice it sometimes does ,
 sometimes doesn ' t . Default is ' default ' , which picks by mode :
 * ' burn ' : False
 * ' do - particles ' : False
 * ' polish ' : True
 dowarn : Bool , optional
 Whether to log a warning if termination results from finishing
 loops rather than from convergence . Default is True .
 Returns
 dictionary
 Dictionary of convergence information . Contains whether the
 optimization has converged ( key ` ` ' converged ' ` ` ) , the values of the
 state after each loop ( key ` ` ' all _ loop _ values ' ` ` ) .
 The values of the state ' s parameters after each part of the
 loop : globals , particles , linemin . If ` ` collect _ stats ` ` is set ,
 then also contains lists of termination dicts from globals ,
 particles , and line minimization ( keys ` ` ' global _ stats ' ` ` ,
 ` ` ' particle _ stats ' ` ` , and ` ` ' line _ stats ` ` ' ) .
 Notes
 Proceeds by alternating between one Levenberg - Marquardt step
 optimizing the globals , one optimizing the particles , and repeating
 until termination .
 In addition , if ` do _ line _ min ` is True , at the end of each loop
 step an additional optimization is tried along the subspaced spanned
 by the steps taken during the last 3 loops . Ideally , this changes the
 convergence from linear to quadratic , but it doesn ' t always do much .
 Each of the 3 options proceed by optimizing as follows :
 * burn : lm . do _ run _ 2 ( ) , lp . do _ run _ 2 ( ) . No psf , 2 loops on lm .
 * do - particles : lp . do _ run _ 2 ( ) , scales for ilm , bkg ' s
 * polish : lm . do _ run _ 2 ( ) , lp . do _ run _ 2 ( ) . Everything , 1 loop each .
 where lm is a globals LMGlobals instance , and lp a
 LMParticleGroupCollection instance ."""

implementation:
# It would be nice if some of these magic # ' s ( region size ,
# num _ eig _ dirs , etc ) were calculated in a good way . FIXME
mode = mode . lower ( )
if mode not in { 'burn' , 'do-particles' , 'polish' } :
    raise ValueError ( 'mode must be one of burn, do-particles, polish' )
# 1 . Setting Defaults
if desc is '' :
    desc = mode + 'ing' if mode != 'do-particles' else 'doing-particles'
eig_update = ( mode != 'do-particles' )
glbl_run_length = 3 if mode == 'do-particles' else 6
glbl_mx_itr = 2 if mode == 'burn' else 1
use_accel = ( mode == 'burn' )
rz_order = int ( rz_order )
if do_line_min == 'default' : # do _ line _ min = ( mode = = ' polish ' )
# temporary fix until we solve the particles - leaving - image issue :
    do_line_min = False
if mode == 'do-particles' : # FIXME explicit params
# We pick some parameters for an overall illumination scale :
    glbl_nms = [ 'ilm-scale' , 'ilm-xy-0-0' , 'bkg-xy-0-0' , 'offset' ]
    # And now , since we have explicit parameters , we check that they
    # are actually in the state :
    glbl_nms = [ nm for nm in glbl_nms if nm in s . params ]
else :
    if mode == 'polish' :
        remove_params = None
    else : # FIXME explicit params
        remove_params = s . get ( 'psf' ) . params
        if ( 'zscale' not in remove_params ) and ( 'zscale' in s . params ) :
            remove_params . append ( 'zscale' )
    glbl_nms = name_globals ( s , remove_params = remove_params )
all_lp_stats = [ ]
all_lm_stats = [ ]
all_line_stats = [ ]
all_loop_values = [ ]
_delta_vals = [ ]
# storing the directions we ' ve moved along for line min
# 2 . Optimize
CLOG . info ( 'Start of loop %d:\t%f' % ( 0 , s . error ) )
for a in range ( n_loop ) :
    start_err = s . error
    start_params = np . copy ( s . state [ s . params ] )
    # 2a . Globals
    # glbl _ dmp = 0.3 if a = = 0 else 3e - 2
    # # # # FIXME we damp degenerate but convenient spaces in the ilm , bkg
    # # # # manually , but we should do it more betterer .
    BAD_DAMP = 1e7
    BAD_LIST = [ [ 'ilm-scale' , BAD_DAMP ] , [ 'ilm-off' , BAD_DAMP ] , [ 'ilm-z-0' , BAD_DAMP ] , [ 'bkg-z-0' , BAD_DAMP ] ]
    glbl_dmp = vectorize_damping ( glbl_nms + [ 'rz' ] * rz_order , damping = 1.0 , increase_list = [ [ 'psf-' , 3e1 ] ] + BAD_LIST )
    if a != 0 or mode != 'do-particles' :
        if partial_log :
            log . set_level ( 'debug' )
        gstats = do_levmarq ( s , glbl_nms , max_iter = glbl_mx_itr , run_length = glbl_run_length , eig_update = eig_update , num_eig_dirs = 10 , eig_update_frequency = 3 , rz_order = rz_order , damping = glbl_dmp , decrease_damp_factor = 10. , use_accel = use_accel , collect_stats = collect_stats , fractol = 0.1 * fractol , max_mem = max_mem )
        if partial_log :
            log . set_level ( 'info' )
        all_lm_stats . append ( gstats )
    if desc is not None :
        states . save ( s , desc = desc )
    CLOG . info ( 'Globals,   loop {}:\t{}' . format ( a , s . error ) )
    all_loop_values . append ( s . values )
    # 2b . Particles
    prtl_dmp = 1.0 if a == 0 else 1e-2
    # For now , I ' m calculating the region size . This might be a bad idea
    # because 1 bad particle can spoil the whole group .
    pstats = do_levmarq_all_particle_groups ( s , region_size = 40 , max_iter = 1 , do_calc_size = True , run_length = 4 , eig_update = False , damping = prtl_dmp , fractol = 0.1 * fractol , collect_stats = collect_stats , max_mem = max_mem , include_rad = include_rad )
    all_lp_stats . append ( pstats )
    if desc is not None :
        states . save ( s , desc = desc )
    CLOG . info ( 'Particles, loop {}:\t{}' . format ( a , s . error ) )
    gc . collect ( )
    all_loop_values . append ( s . values )
    # 2c . Line min ?
    end_params = np . copy ( s . state [ s . params ] )
    _delta_vals . append ( start_params - end_params )
    if do_line_min :
        all_line_stats . append ( do_levmarq_n_directions ( s , _delta_vals [ - 3 : ] , collect_stats = collect_stats ) )
        if desc is not None :
            states . save ( s , desc = desc )
        CLOG . info ( 'Line min., loop {}:\t{}' . format ( a , s . error ) )
        all_loop_values . append ( s . values )
    # 2d . terminate ?
    new_err = s . error
    derr = start_err - new_err
    dobreak = ( derr / new_err < fractol ) or ( derr < errtol )
    if dobreak :
        break
if dowarn and ( not dobreak ) :
    CLOG . warn ( 'burn() did not converge; consider re-running' )
d = { 'converged' : dobreak , 'all_loop_values' : all_loop_values }
if collect_stats :
    d . update ( { 'global_stats' : all_lm_stats , 'particle_stats' : all_lp_stats , 'line_stats' : all_line_stats } )
return d


signature:
def finish ( s , desc = 'finish' , n_loop = 4 , max_mem = 1e9 , separate_psf = True , fractol = 1e-7 , errtol = 1e-3 , dowarn = True ) :
    """Crawls slowly to the minimum - cost state .
 Blocks the global parameters into small enough sections such that each
 can be optimized separately while including all the pixels ( i . e . no
 decimation ) . Optimizes the globals , then the psf separately if desired ,
 then particles , then a line minimization along the step direction to
 speed up convergence .
 Parameters
 s : : class : ` peri . states . ImageState `
 The state to optimize
 desc : string , optional
 Description to append to the states . save ( ) call every loop .
 Set to ` None ` to avoid saving . Default is ` ' finish ' ` .
 n _ loop : Int , optional
 The number of times to loop over in the optimizer . Default is 4.
 max _ mem : Numeric , optional
 The maximum amount of memory allowed for the optimizers ' J ' s ,
 for both particles & globals . Default is 1e9.
 separate _ psf : Bool , optional
 If True , does the psf optimization separately from the rest of
 the globals , since the psf has a more tortuous fit landscape .
 Default is True .
 fractol : Float , optional
 Fractional change in error at which to terminate . Default 1e - 4
 errtol : Float , optional
 Absolute change in error at which to terminate . Default 1e - 2
 dowarn : Bool , optional
 Whether to log a warning if termination results from finishing
 loops rather than from convergence . Default is True .
 Returns
 dictionary
 Information about the optimization . Has two keys : ` ` ' converged ' ` ` ,
 a Bool which of whether optimization stopped due to convergence
 ( True ) or due to max number of iterations ( False ) , and
 ` ` ' loop _ values ' ` ` , a [ n _ loop + 1 , N ] ` ` numpy . ndarray ` ` of the
 state ' s values , at the start of optimization and at the end of
 each loop , before the line minimization ."""

implementation:
values = [ np . copy ( s . state [ s . params ] ) ]
remove_params = s . get ( 'psf' ) . params if separate_psf else None
# FIXME explicit params
global_params = name_globals ( s , remove_params = remove_params )
# FIXME this could be done much better , since much of the globals such
# as the ilm are local . Could be done with sparse matrices and / or taking
# nearby globals in a group and using the update tile only as the slicer ,
# rather than the full residuals .
gs = np . floor ( max_mem / s . residuals . nbytes ) . astype ( 'int' )
groups = [ global_params [ a : a + gs ] for a in range ( 0 , len ( global_params ) , gs ) ]
CLOG . info ( 'Start  ``finish``:\t{}' . format ( s . error ) )
for a in range ( n_loop ) :
    start_err = s . error
    # 1 . Min globals :
    for g in groups :
        do_levmarq ( s , g , damping = 0.1 , decrease_damp_factor = 20. , max_iter = 1 , max_mem = max_mem , eig_update = False )
    if separate_psf :
        do_levmarq ( s , remove_params , max_mem = max_mem , max_iter = 4 , eig_update = False )
    CLOG . info ( 'Globals,   loop {}:\t{}' . format ( a , s . error ) )
    if desc is not None :
        states . save ( s , desc = desc )
    # 2 . Min particles
    do_levmarq_all_particle_groups ( s , max_iter = 1 , max_mem = max_mem )
    CLOG . info ( 'Particles, loop {}:\t{}' . format ( a , s . error ) )
    if desc is not None :
        states . save ( s , desc = desc )
    # 3 . Append vals , line min :
    values . append ( np . copy ( s . state [ s . params ] ) )
    # dv = ( np . array ( values [ 1 : ] ) - np . array ( values [ 0 ] ) ) [ - 3 : ]
    # do _ levmarq _ n _ directions ( s , dv , damping = 1e - 2 , max _ iter = 2 , errtol = 3e - 4)
    # CLOG . info ( ' Line min . , loop { } : \ t { } ' . format ( a , s . error ) )
    # if desc is not None :
    # states . save ( s , desc = desc )
    # 4 . terminate ?
    new_err = s . error
    derr = start_err - new_err
    dobreak = ( derr / new_err < fractol ) or ( derr < errtol )
    if dobreak :
        break
if dowarn and ( not dobreak ) :
    CLOG . warn ( 'finish() did not converge; consider re-running' )
return { 'converged' : dobreak , 'loop_values' : np . array ( values ) }


signature:
def fit_comp ( new_comp , old_comp , ** kwargs ) :
    """Fits a new component to an old component
 Calls do _ levmarq to match the . get ( ) fields of the two objects . The
 parameters of new _ comp are modified in place .
 Parameters
 new _ comp : : class : ` peri . comps . comp `
 The new object , whose parameters to update to fit the field of
 ` old _ comp ` . Must have a . get ( ) attribute which returns an ndarray
 old _ comp : peri . comp
 The old ilm to match to .
 Other Parameters
 Any keyword arguments to be passed to the optimizer LMGlobals
 through do _ levmarq .
 See Also
 do _ levmarq : Levenberg - Marquardt minimization using a random subset
 of the image pixels ."""

implementation:
# resetting the category to ilm :
new_cat = new_comp . category
new_comp . category = 'ilm'
fake_s = states . ImageState ( Image ( old_comp . get ( ) . copy ( ) ) , [ new_comp ] , pad = 0 , mdl = mdl . SmoothFieldModel ( ) )
do_levmarq ( fake_s , new_comp . params , ** kwargs )
new_comp . category = new_cat


signature:
def reset ( self , new_damping = None ) :
    """Keeps all user supplied options the same , but resets counters etc ."""

implementation:
self . _num_iter = 0
self . _inner_run_counter = 0
self . _J_update_counter = self . update_J_frequency
self . _fresh_JTJ = False
self . _has_run = False
if new_damping is not None :
    self . damping = np . array ( new_damping ) . astype ( 'float' )
self . _set_err_paramvals ( )


signature:
def do_run_1 ( self ) :
    """LM run , evaluating 1 step at a time .
 Broyden or eigendirection updates replace full - J updates until
 a full - J update occurs . Does not run with the calculated J ( no
 internal run ) ."""

implementation:
while not self . check_terminate ( ) :
    self . _has_run = True
    self . _run1 ( )
    self . _num_iter += 1 ;
    self . _inner_run_counter += 1


signature:
def _run1 ( self ) :
    """workhorse for do _ run _ 1"""

implementation:
if self . check_update_J ( ) :
    self . update_J ( )
else :
    if self . check_Broyden_J ( ) :
        self . update_Broyden_J ( )
    if self . check_update_eig_J ( ) :
        self . update_eig_J ( )
# 1 . Assuming that J starts updated :
delta_vals = self . find_LM_updates ( self . calc_grad ( ) )
# 2 . Increase damping until we get a good step :
er1 = self . update_function ( self . param_vals + delta_vals )
good_step = ( find_best_step ( [ self . error , er1 ] ) == 1 )
if not good_step :
    er0 = self . update_function ( self . param_vals )
    if np . abs ( er0 - self . error ) / er0 > 1e-7 :
        raise RuntimeError ( 'Function updates are not exact.' )
    CLOG . debug ( 'Bad step, increasing damping' )
    CLOG . debug ( '\t\t%f\t%f' % ( self . error , er1 ) )
    grad = self . calc_grad ( )
    for _try in range ( self . _max_inner_loop ) :
        self . increase_damping ( )
        delta_vals = self . find_LM_updates ( grad )
        er1 = self . update_function ( self . param_vals + delta_vals )
        good_step = ( find_best_step ( [ self . error , er1 ] ) == 1 )
        if good_step :
            break
    else :
        er0 = self . update_function ( self . param_vals )
        CLOG . warn ( 'Stuck!' )
        if np . abs ( er0 - self . error ) / er0 > 1e-7 :
            raise RuntimeError ( 'Function updates are not exact.' )
# state is updated , now params :
if good_step :
    self . _last_error = self . error
    self . error = er1
    CLOG . debug ( 'Good step\t%f\t%f' % ( self . _last_error , self . error ) )
    self . update_param_vals ( delta_vals , incremental = True )
    self . decrease_damping ( )


signature:
def do_run_2 ( self ) :
    """LM run evaluating 2 steps ( damped and not ) and choosing the best .
 After finding the best of 2 steps , runs with that damping + Broyden
 or eigendirection updates , until deciding to do a full - J update .
 Only changes damping after full - J updates ."""

implementation:
while not self . check_terminate ( ) :
    self . _has_run = True
    self . _run2 ( )
    self . _num_iter += 1


signature:
def _run2 ( self ) :
    """Workhorse for do _ run _ 2"""

implementation:
if self . check_update_J ( ) :
    self . update_J ( )
else :
    if self . check_Broyden_J ( ) :
        self . update_Broyden_J ( )
    if self . check_update_eig_J ( ) :
        self . update_eig_J ( )
# 0 . Find _ last _ residuals , _ last _ error , etc :
_last_residuals = self . calc_residuals ( ) . copy ( )
_last_error = 1 * self . error
_last_vals = self . param_vals . copy ( )
# 1 . Calculate 2 possible steps
delta_params_1 = self . find_LM_updates ( self . calc_grad ( ) , do_correct_damping = False )
self . decrease_damping ( )
delta_params_2 = self . find_LM_updates ( self . calc_grad ( ) , do_correct_damping = False )
self . decrease_damping ( undo_decrease = True )
# 2 . Check which step is best :
er1 = self . update_function ( self . param_vals + delta_params_1 )
er2 = self . update_function ( self . param_vals + delta_params_2 )
triplet = ( self . error , er1 , er2 )
best_step = find_best_step ( triplet )
if best_step == 0 : # Both bad steps , put back & increase damping :
    _ = self . update_function ( self . param_vals . copy ( ) )
    grad = self . calc_grad ( )
    CLOG . debug ( 'Bad step, increasing damping' )
    CLOG . debug ( '%f\t%f\t%f' % triplet )
    for _try in range ( self . _max_inner_loop ) :
        self . increase_damping ( )
        delta_vals = self . find_LM_updates ( grad )
        er_new = self . update_function ( self . param_vals + delta_vals )
        good_step = er_new < self . error
        if good_step : # Update params , error , break :
            self . update_param_vals ( delta_vals , incremental = True )
            self . error = er_new
            CLOG . debug ( 'Sufficiently increased damping' )
            CLOG . debug ( '%f\t%f' % ( triplet [ 0 ] , self . error ) )
            break
    else : # for - break - else
    # Throw a warning , put back the parameters
        CLOG . warn ( 'Stuck!' )
        self . error = self . update_function ( self . param_vals . copy ( ) )
elif best_step == 1 : # er1 < = er2:
    good_step = True
    CLOG . debug ( 'Good step, same damping' )
    CLOG . debug ( '%f\t%f\t%f' % triplet )
    # Update to er1 params :
    er1_1 = self . update_function ( self . param_vals + delta_params_1 )
    if np . abs ( er1_1 - er1 ) > 1e-6 :
        raise RuntimeError ( 'Function updates are not exact.' )
    self . update_param_vals ( delta_params_1 , incremental = True )
    self . error = er1
elif best_step == 2 : # er2 < er1:
    good_step = True
    self . error = er2
    CLOG . debug ( 'Good step, decreasing damping' )
    CLOG . debug ( '%f\t%f\t%f' % triplet )
    # - we ' re already at the correct parameters
    self . update_param_vals ( delta_params_2 , incremental = True )
    self . decrease_damping ( )
# 3 . Run with current J , damping ; update what we need to : :
if good_step :
    self . _last_residuals = _last_residuals
    self . _last_error = _last_error
    self . _last_vals = _last_vals
    self . error
    self . do_internal_run ( initial_count = 1 )


signature:
def do_internal_run ( self , initial_count = 0 , subblock = None , update_derr = True ) :
    """Takes more steps without calculating J again .
 Given a fixed damping , J , JTJ , iterates calculating steps , with
 optional Broyden or eigendirection updates . Iterates either until
 a bad step is taken or for self . run _ length times .
 Called internally by do _ run _ 2 ( ) but is also useful on its own .
 Parameters
 initial _ count : Int , optional
 The initial count of the run . Default is 0 . Increasing from
 0 effectively temporarily decreases run _ length .
 subblock : None or np . ndarray of bools , optional
 If not None , a boolean mask which determines which sub -
 block of parameters to run over . Default is None , i . e .
 all the parameters .
 update _ derr : Bool , optional
 Set to False to not update the variable that determines
 delta _ err , preventing premature termination through errtol .
 Notes
 It might be good to do something similar to update _ derr with the
 parameter values , but this is trickier because of Broyden updates
 and _ fresh _ J ."""

implementation:
self . _inner_run_counter = initial_count ;
good_step = True
n_good_steps = 0
CLOG . debug ( 'Running...' )
_last_residuals = self . calc_residuals ( ) . copy ( )
while ( ( self . _inner_run_counter < self . run_length ) & good_step & ( not self . check_terminate ( ) ) ) : # 1 . Checking if we update J
    if self . check_Broyden_J ( ) and self . _inner_run_counter != 0 :
        self . update_Broyden_J ( )
    if self . check_update_eig_J ( ) and self . _inner_run_counter != 0 :
        self . update_eig_J ( )
    # 2 . Getting parameters , error
    er0 = 1 * self . error
    delta_vals = self . find_LM_updates ( self . calc_grad ( ) , do_correct_damping = False , subblock = subblock )
    er1 = self . update_function ( self . param_vals + delta_vals )
    good_step = er1 < er0
    if good_step :
        n_good_steps += 1
        CLOG . debug ( '%f\t%f' % ( er0 , er1 ) )
        # Updating :
        self . update_param_vals ( delta_vals , incremental = True )
        self . _last_residuals = _last_residuals . copy ( )
        if update_derr :
            self . _last_error = er0
        self . error = er1
        _last_residuals = self . calc_residuals ( ) . copy ( )
    else :
        er0_0 = self . update_function ( self . param_vals )
        CLOG . debug ( 'Bad step!' )
        if np . abs ( er0 - er0_0 ) > 1e-6 :
            raise RuntimeError ( 'Function updates are not exact.' )
    self . _inner_run_counter += 1
return n_good_steps


signature:
def find_LM_updates ( self , grad , do_correct_damping = True , subblock = None ) :
    """Calculates LM updates , with or without the acceleration correction .
 Parameters
 grad : numpy . ndarray
 The gradient of the model cost .
 do _ correct _ damping : Bool , optional
 If ` self . use _ accel ` , then set to True to correct damping
 if the acceleration correction is too big . Default is True
 Does nothing is ` self . use _ accel ` is False
 subblock : slice , numpy . ndarray , or None , optional
 Set to a slice or a valide numpy . ndarray to use only a
 certain subset of the parameters . Default is None , i . e .
 use all the parameters .
 Returns
 delta : numpy . ndarray
 The Levenberg - Marquadt step , relative to the old
 parameters . Size is always self . param _ vals . size ."""

implementation:
if subblock is not None :
    if ( subblock . sum ( ) == 0 ) or ( subblock . size == 0 ) :
        CLOG . fatal ( 'Empty subblock in find_LM_updates' )
        raise ValueError ( 'Empty sub-block' )
    j = self . J [ subblock ]
    JTJ = np . dot ( j , j . T )
    damped_JTJ = self . _calc_damped_jtj ( JTJ , subblock = subblock )
    grad = grad [ subblock ]
    # select the subblock of the grad
else :
    damped_JTJ = self . _calc_damped_jtj ( self . JTJ , subblock = subblock )
delta = self . _calc_lm_step ( damped_JTJ , grad , subblock = subblock )
if self . use_accel :
    accel_correction = self . calc_accel_correction ( damped_JTJ , delta )
    nrm_d0 = np . sqrt ( np . sum ( delta ** 2 ) )
    nrm_corr = np . sqrt ( np . sum ( accel_correction ** 2 ) )
    CLOG . debug ( '|correction| / |LM step|\t%e' % ( nrm_corr / nrm_d0 ) )
    if nrm_corr / nrm_d0 < self . max_accel_correction :
        delta += accel_correction
    elif do_correct_damping :
        CLOG . debug ( 'Untrustworthy step! Increasing damping...' )
        self . increase_damping ( )
        damped_JTJ = self . _calc_damped_jtj ( self . JTJ , subblock = subblock )
        delta = self . _calc_lm_step ( damped_JTJ , grad , subblock = subblock )
if np . any ( np . isnan ( delta ) ) :
    CLOG . fatal ( 'Calculated steps have nans!?' )
    raise FloatingPointError ( 'Calculated steps have nans!?' )
return delta


signature:
def _calc_lm_step ( self , damped_JTJ , grad , subblock = None ) :
    """Calculates a Levenberg - Marquard step w / o acceleration"""

implementation:
delta0 , res , rank , s = np . linalg . lstsq ( damped_JTJ , - 0.5 * grad , rcond = self . min_eigval )
if self . _fresh_JTJ :
    CLOG . debug ( '%d degenerate of %d total directions' % ( delta0 . size - rank , delta0 . size ) )
if subblock is not None :
    delta = np . zeros ( self . J . shape [ 0 ] )
    delta [ subblock ] = delta0
else :
    delta = delta0 . copy ( )
return delta


signature:
def update_param_vals ( self , new_vals , incremental = False ) :
    """Updates the current set of parameter values and previous values ,
 sets a flag to re - calculate J .
 Parameters
 new _ vals : numpy . ndarray
 The new values to update to
 incremental : Bool , optional
 Set to True to make it an incremental update relative
 to the old parameters . Default is False"""

implementation:
self . _last_vals = self . param_vals . copy ( )
if incremental :
    self . param_vals += new_vals
else :
    self . param_vals = new_vals . copy ( )
# And we ' ve updated , so JTJ is no longer valid :
self . _fresh_JTJ = False


signature:
def find_expected_error ( self , delta_params = 'calc' ) :
    """Returns the error expected after an update if the model were linear .
 Parameters
 delta _ params : { numpy . ndarray , ' calc ' , or ' perfect ' } , optional
 The relative change in parameters . If ' calc ' , uses update
 calculated from the current damping , J , etc ; if ' perfect ' ,
 uses the update calculated with zero damping .
 Returns
 numpy . float64
 The expected error after the update with ` delta _ params `"""

implementation:
grad = self . calc_grad ( )
if list ( delta_params ) in [ list ( 'calc' ) , list ( 'perfect' ) ] :
    jtj = ( self . JTJ if delta_params == 'perfect' else self . _calc_damped_jtj ( self . JTJ ) )
    delta_params = self . _calc_lm_step ( jtj , self . calc_grad ( ) )
# If the model were linear , then the cost would be quadratic ,
# with Hessian 2 * ` self . JTJ ` and gradient ` grad `
expected_error = ( self . error + np . dot ( grad , delta_params ) + np . dot ( np . dot ( self . JTJ , delta_params ) , delta_params ) )
return expected_error


signature:
def calc_model_cosine ( self , decimate = None , mode = 'err' ) :
    """Calculates the cosine of the residuals with the model .
 Parameters
 decimate : Int or None , optional
 Decimate the residuals by ` decimate ` pixels . If None , no
 decimation is used . Valid only with mode = ' svd ' . Default
 is None
 mode : { ' svd ' , ' err ' }
 Which mode to use ; see Notes section . Default is ' err ' .
 Returns
 abs _ cos : numpy . float64
 The absolute value of the model cosine .
 Notes
 The model cosine is defined in terms of the geometric view of
 curve - fitting , as a model manifold embedded in a high - dimensional
 space . The model cosine is the cosine of the residuals vector
 with its projection on the tangent space : : math : ` cos ( phi ) = | P ^ T r | / | r | `
 where : math : ` P ^ T ` is the projection operator onto the model manifold
 and : math : ` r ` the residuals . This can be calculated two ways : By
 calculating the projection operator P directly with SVD ( mode = ` svd ` ) ,
 or by using the expected error if the model were linear to calculate
 a model sine first ( mode = ` err ` ) . Since the SVD of a large matrix is
 slow , mode = ` err ` is faster .
 ` decimate ` allows for every nth pixel only to be counted in the
 SVD matrix of J for speed . While this is n x faster , it is
 considerably less accurate , so the default is no decimation ."""

implementation:
if mode == 'svd' :
    slicer = slice ( 0 , None , decimate )
    # 1 . Calculate projection term
    u , sig , v = np . linalg . svd ( self . J [ : , slicer ] , full_matrices = False )
    # slow part
    # p = np . dot ( v . T , v ) - memory error , so term - by - term
    r = self . calc_residuals ( ) [ slicer ]
    abs_r = np . sqrt ( ( r * r ) . sum ( ) )
    v_r = np . dot ( v , r / abs_r )
    projected = np . dot ( v . T , v_r )
    abs_cos = np . sqrt ( ( projected * projected ) . sum ( ) )
elif mode == 'err' :
    expected_error = self . find_expected_error ( delta_params = 'perfect' )
    model_sine_2 = expected_error / self . error
    # error = distance ^ 2
    abs_cos = np . sqrt ( 1 - model_sine_2 )
else :
    raise ValueError ( 'mode must be one of `svd`, `err`' )
return abs_cos


signature:
def get_termination_stats ( self , get_cos = True ) :
    """Returns a dict of termination statistics
 Parameters
 get _ cos : Bool , optional
 Whether or not to calcualte the cosine of the residuals
 with the tangent plane of the model using the current J .
 The calculation may take some time . Default is True
 Returns
 dict
 Has keys
 delta _ vals : The last change in parameter values .
 delta _ err : The last change in the error .
 exp _ err : The expected ( last ) change in the error .
 frac _ err : The fractional change in the error .
 num _ iter : The number of iterations completed .
 error : The current error ."""

implementation:
delta_vals = self . _last_vals - self . param_vals
delta_err = self . _last_error - self . error
frac_err = delta_err / self . error
to_return = { 'delta_vals' : delta_vals , 'delta_err' : delta_err , 'num_iter' : 1 * self . _num_iter , 'frac_err' : frac_err , 'error' : self . error , 'exp_err' : self . _exp_err }
if get_cos :
    model_cosine = self . calc_model_cosine ( )
    to_return . update ( { 'model_cosine' : model_cosine } )
return to_return


signature:
def check_completion ( self ) :
    """Returns a Bool of whether the algorithm has found a satisfactory minimum"""

implementation:
terminate = False
term_dict = self . get_termination_stats ( get_cos = self . costol is not None )
terminate |= np . all ( np . abs ( term_dict [ 'delta_vals' ] ) < self . paramtol )
terminate |= ( term_dict [ 'delta_err' ] < self . errtol )
terminate |= ( term_dict [ 'exp_err' ] < self . exptol )
terminate |= ( term_dict [ 'frac_err' ] < self . fractol )
if self . costol is not None :
    terminate |= ( curcos < term_dict [ 'model_cosine' ] )
return terminate


signature:
def check_terminate ( self ) :
    """Returns a Bool of whether to terminate .
 Checks whether a satisfactory minimum has been found or whether
 too many iterations have occurred ."""

implementation:
if not self . _has_run :
    return False
else : # 1-3 . errtol , paramtol , model cosine low enough ?
    terminate = self . check_completion ( )
    # 4 . too many iterations ? ?
    terminate |= ( self . _num_iter >= self . max_iter )
    return terminate


signature:
def check_update_J ( self ) :
    """Checks if the full J should be updated .
 Right now , just updates after update _ J _ frequency loops"""

implementation:
self . _J_update_counter += 1
update = self . _J_update_counter >= self . update_J_frequency
return update & ( not self . _fresh_JTJ )


signature:
def update_J ( self ) :
    """Updates J , JTJ , and internal counters ."""

implementation:
self . calc_J ( )
# np . dot ( j , j . T ) is slightly faster but 2x as much mem
step = np . ceil ( 1e-2 * self . J . shape [ 1 ] ) . astype ( 'int' )
# 1 % more mem . . .
self . JTJ = low_mem_sq ( self . J , step = step )
# copies still , since J is not C - ordered but a slice of j _ e . . .
# doing self . J . copy ( ) works but takes 2x as much ram . .
self . _fresh_JTJ = True
self . _J_update_counter = 0
if np . any ( np . isnan ( self . JTJ ) ) :
    raise FloatingPointError ( 'J, JTJ have nans.' )
# Update self . _ exp _ err
self . _exp_err = self . error - self . find_expected_error ( delta_params = 'perfect' )


signature:
def calc_grad ( self ) :
    """The gradient of the cost w . r . t . the parameters ."""

implementation:
residuals = self . calc_residuals ( )
return 2 * np . dot ( self . J , residuals )


signature:
def _rank_1_J_update ( self , direction , values ) :
    """Does J + = np . outer ( direction , new _ values - old _ values ) without
 using lots of memory"""

implementation:
vals_to_sub = np . dot ( direction , self . J )
delta_vals = values - vals_to_sub
for a in range ( direction . size ) :
    self . J [ a ] += direction [ a ] * delta_vals


signature:
def update_Broyden_J ( self ) :
    """Execute a Broyden update of J"""

implementation:
CLOG . debug ( 'Broyden update.' )
delta_vals = self . param_vals - self . _last_vals
delta_residuals = self . calc_residuals ( ) - self . _last_residuals
nrm = np . sqrt ( np . dot ( delta_vals , delta_vals ) )
direction = delta_vals / nrm
vals = delta_residuals / nrm
self . _rank_1_J_update ( direction , vals )
self . JTJ = np . dot ( self . J , self . J . T )


signature:
def update_eig_J ( self ) :
    """Execute an eigen update of J"""

implementation:
CLOG . debug ( 'Eigen update.' )
vls , vcs = np . linalg . eigh ( self . JTJ )
res0 = self . calc_residuals ( )
for a in range ( min ( [ self . num_eig_dirs , vls . size ] ) ) : # 1 . Finding stiff directions
    stif_dir = vcs [ - ( a + 1 ) ]
    # already normalized
    # 2 . Evaluating derivative along that direction , we ' ll use dl = 5e - 4:
    dl = self . eig_dl
    # 1e - 5
    _ = self . update_function ( self . param_vals + dl * stif_dir )
    res1 = self . calc_residuals ( )
    # 3 . Updating
    grad_stif = ( res1 - res0 ) / dl
    self . _rank_1_J_update ( stif_dir , grad_stif )
self . JTJ = np . dot ( self . J , self . J . T )
# Putting the parameters back :
_ = self . update_function ( self . param_vals )


signature:
def calc_accel_correction ( self , damped_JTJ , delta0 ) :
    """Geodesic acceleration correction to the LM step .
 Parameters
 damped _ JTJ : numpy . ndarray
 The damped JTJ used to calculate the initial step .
 delta0 : numpy . ndarray
 The initial LM step .
 Returns
 corr : numpy . ndarray
 The correction to the original LM step ."""

implementation:
# Get the derivative :
_ = self . update_function ( self . param_vals )
rm0 = self . calc_residuals ( ) . copy ( )
_ = self . update_function ( self . param_vals + delta0 )
rm1 = self . calc_residuals ( ) . copy ( )
_ = self . update_function ( self . param_vals - delta0 )
rm2 = self . calc_residuals ( ) . copy ( )
der2 = ( rm2 + rm1 - 2 * rm0 )
corr , res , rank , s = np . linalg . lstsq ( damped_JTJ , np . dot ( self . J , der2 ) , rcond = self . min_eigval )
corr *= - 0.5
return corr


signature:
def update_select_J ( self , blk ) :
    """Updates J only for certain parameters , described by the boolean
 mask ` blk ` ."""

implementation:
p0 = self . param_vals . copy ( )
self . update_function ( p0 )
# in case things are not put back . . .
r0 = self . calc_residuals ( ) . copy ( )
dl = np . zeros ( p0 . size , dtype = 'float' )
blk_J = [ ]
for i in np . nonzero ( blk ) [ 0 ] :
    dl *= 0 ;
    dl [ i ] = self . eig_dl
    self . update_function ( p0 + dl )
    r1 = self . calc_residuals ( ) . copy ( )
    blk_J . append ( ( r1 - r0 ) / self . eig_dl )
self . J [ blk ] = np . array ( blk_J )
self . update_function ( p0 )
# Then we also need to update JTJ :
self . JTJ = np . dot ( self . J , self . J . T )
if np . any ( np . isnan ( self . J ) ) or np . any ( np . isnan ( self . JTJ ) ) :
    raise FloatingPointError ( 'J, JTJ have nans.' )


signature:
def _set_err_paramvals ( self ) :
    """Must update :
 self . error , self . _ last _ error , self . param _ vals , self . _ last _ vals"""

implementation:
# self . param _ vals = p0 # sloppy . . .
self . _last_vals = self . param_vals . copy ( )
self . error = self . update_function ( self . param_vals )
self . _last_error = ( 1 + 2 * self . fractol ) * self . error


signature:
def calc_J ( self ) :
    """Updates self . J , returns nothing"""

implementation:
del self . J
self . J = np . zeros ( [ self . param_vals . size , self . data . size ] )
dp = np . zeros_like ( self . param_vals )
f0 = self . model . copy ( )
for a in range ( self . param_vals . size ) :
    dp *= 0
    dp [ a ] = self . dl [ a ]
    f1 = self . func ( self . param_vals + dp , * self . func_args , ** self . func_kwargs )
    grad_func = ( f1 - f0 ) / dp [ a ]
    # J = grad ( residuals ) = - grad ( model )
    self . J [ a ] = - grad_func


signature:
def update_function ( self , param_vals ) :
    """Takes an array param _ vals , updates function , returns the new error"""

implementation:
self . model = self . func ( param_vals , * self . func_args , ** self . func_kwargs )
d = self . calc_residuals ( )
return np . dot ( d . flat , d . flat )


signature:
def update_function ( self , param_vals ) :
    """Updates the opt _ obj , returns new error ."""

implementation:
self . opt_obj . update_function ( param_vals )
return self . opt_obj . get_error ( )


signature:
def update_function ( self , param_vals ) :
    """Updates with param _ vals [ i ] = distance from self . p0 along self . direction [ i ] ."""

implementation:
dp = np . zeros ( self . p0 . size )
for a in range ( param_vals . size ) :
    dp += param_vals [ a ] * self . directions [ a ]
self . state . update ( self . state . params , self . p0 + dp )
self . param_vals [ : ] = param_vals
return None


signature:
def calc_J ( self ) :
    """Calculates J along the direction ."""

implementation:
r0 = self . state . residuals . copy ( ) . ravel ( )
dl = np . zeros ( self . param_vals . size )
p0 = self . param_vals . copy ( )
J = [ ]
for a in range ( self . param_vals . size ) :
    dl *= 0
    dl [ a ] += self . dl
    self . update_function ( p0 + dl )
    r1 = self . state . residuals . copy ( ) . ravel ( )
    J . append ( ( r1 - r0 ) / self . dl )
self . update_function ( p0 )
return np . array ( J )


signature:
def update_select_J ( self , blk ) :
    """Updates J only for certain parameters , described by the boolean
 mask blk ."""

implementation:
self . update_function ( self . param_vals )
params = np . array ( self . param_names ) [ blk ] . tolist ( )
blk_J = - self . state . gradmodel ( params = params , inds = self . _inds , flat = False )
self . J [ blk ] = blk_J
# Then we also need to update JTJ :
self . JTJ = np . dot ( self . J , self . J . T )
if np . any ( np . isnan ( self . J ) ) or np . any ( np . isnan ( self . JTJ ) ) :
    raise FloatingPointError ( 'J, JTJ have nans.' )


signature:
def find_expected_error ( self , delta_params = 'calc' , adjust = True ) :
    """Returns the error expected after an update if the model were linear .
 Parameters
 delta _ params : { numpy . ndarray , ' calc ' , or ' perfect ' } , optional
 The relative change in parameters . If ' calc ' , uses update
 calculated from the current damping , J , etc ; if ' perfect ' ,
 uses the update calculated with zero damping .
 Returns
 numpy . float64
 The expected error after the update with ` delta _ params `"""

implementation:
expected_error = super ( LMGlobals , self ) . find_expected_error ( delta_params = delta_params )
if adjust : # adjust for num _ pix
    derr = ( expected_error - self . error ) * ( self . state . residuals . size / float ( self . num_pix ) )
    expected_error = self . error + derr
return expected_error


signature:
def calc_model_cosine ( self , decimate = None , mode = 'err' ) :
    """Calculates the cosine of the residuals with the model .
 Parameters
 decimate : Int or None , optional
 Decimate the residuals by ` decimate ` pixels . If None , no
 decimation is used . Valid only with mode = ' svd ' . Default
 is None
 mode : { ' svd ' , ' err ' }
 Which mode to use ; see Notes section . Default is ' err ' .
 Returns
 abs _ cos : numpy . float64
 The absolute value of the model cosine .
 Notes
 The model cosine is defined in terms of the geometric view of
 curve - fitting , as a model manifold embedded in a high - dimensional
 space . The model cosine is the cosine of the residuals vector
 with its projection on the tangent space : : math : ` cos ( phi ) = | P ^ T r | / | r | `
 where : math : ` P ^ T ` is the projection operator onto the model manifold
 and : math : ` r ` the residuals . This can be calculated two ways : By
 calculating the projection operator P directly with SVD ( mode = ` svd ` ) ,
 or by using the expected error if the model were linear to calculate
 a model sine first ( mode = ` err ` ) . Since the SVD of a large matrix is
 slow , mode = ` err ` is faster .
 ` decimate ` allows for every nth pixel only to be counted in the
 SVD matrix of J for speed . While this is n x faster , it is
 considerably less accurate , so the default is no decimation ."""

implementation:
# we calculate the model cosine only in the data space of the
# sampled indices
if mode == 'err' :
    expected_error = self . find_expected_error ( delta_params = 'perfect' , adjust = False )
    derr = self . error - expected_error
    residuals_err = lambda r : np . dot ( r , r ) . sum ( )
    current_partial_error = residuals_err ( self . calc_residuals ( ) )
    expected_partial_error = current_partial_error - derr
    model_sine_2 = expected_partial_error / current_partial_error
    abs_cos = np . sqrt ( 1 - model_sine_2 )
else : # superclass is fine
    abs_cos = super ( self . __class__ , self ) . calc_model_cosine ( decimate = decimate , mode = mode )
return abs_cos


implementation:
if self . _fresh_JTJ :
    return self . _graderr
else :
    residuals = self . calc_residuals ( )
    return 2 * np . dot ( self . J , residuals )


signature:
def reset ( self , new_region_size = None , do_calc_size = True , new_damping = None , new_max_mem = None ) :
    """Resets the particle groups and optionally the region size and damping .
 Parameters
 new _ region _ size : : Int or 3 - element list - like of ints , optional
 The region size for sub - blocking particles . Default is 40
 do _ calc _ size : Bool , optional
 If True , calculates the region size internally based on
 the maximum allowed memory . Default is True
 new _ damping : Float or None , optional
 The new damping of the optimizer . Set to None to leave
 as the default for LMParticles . Default is None .
 new _ max _ mem : Numeric , optional
 The maximum allowed memory for J to occupy . Default is 1e9"""

implementation:
if new_region_size is not None :
    self . region_size = new_region_size
if new_max_mem != None :
    self . max_mem = new_max_mem
if do_calc_size :
    self . region_size = calc_particle_group_region_size ( self . state , region_size = self . region_size , max_mem = self . max_mem )
self . stats = [ ]
self . particle_groups = separate_particles_into_groups ( self . state , self . region_size , doshift = 'rand' )
if new_damping is not None :
    self . _kwargs . update ( { 'damping' : new_damping } )
if self . save_J :
    if len ( self . particle_groups ) > 90 :
        CLOG . warn ( 'Attempting to create many open files. Consider increasing max_mem and/or region_size to avoid crashes.' )
    self . _tempfiles = [ ]
    self . _has_saved_J = [ ]
    for a in range ( len ( self . particle_groups ) ) : # TemporaryFile is automatically deleted
        for _ in [ 'j' , 'tile' ] :
            self . _tempfiles . append ( tempfile . TemporaryFile ( dir = os . getcwd ( ) ) )
        self . _has_saved_J . append ( False )


signature:
def _do_run ( self , mode = '1' ) :
    """workhorse for the self . do _ run _ xx methods ."""

implementation:
for a in range ( len ( self . particle_groups ) ) :
    group = self . particle_groups [ a ]
    lp = LMParticles ( self . state , group , ** self . _kwargs )
    if mode == 'internal' :
        lp . J , lp . JTJ , lp . _dif_tile = self . _load_j_diftile ( a )
    if mode == '1' :
        lp . do_run_1 ( )
    if mode == '2' :
        lp . do_run_2 ( )
    if mode == 'internal' :
        lp . do_internal_run ( )
    self . stats . append ( lp . get_termination_stats ( get_cos = self . get_cos ) )
    if self . save_J and ( mode != 'internal' ) :
        self . _dump_j_diftile ( a , lp . J , lp . _dif_tile )
        self . _has_saved_J [ a ] = True


signature:
def do_internal_run ( self ) :
    """Calls LMParticles . do _ internal _ run for each group of particles ."""

implementation:
if not self . save_J :
    raise RuntimeError ( 'self.save_J=True required for do_internal_run()' )
if not np . all ( self . _has_saved_J ) :
    raise RuntimeError ( 'J, JTJ have not been pre-computed. Call do_run_1 or do_run_2' )
self . _do_run ( mode = 'internal' )


signature:
def reset ( self ) :
    """Resets the initial radii used for updating the particles . Call
 if any of the particle radii or positions have been changed
 external to the augmented state ."""

implementation:
inds = list ( range ( self . state . obj_get_positions ( ) . shape [ 0 ] ) )
self . _rad_nms = self . state . param_particle_rad ( inds )
self . _pos_nms = self . state . param_particle_pos ( inds )
self . _initial_rad = np . copy ( self . state . state [ self . _rad_nms ] )
self . _initial_pos = np . copy ( self . state . state [ self . _pos_nms ] ) . reshape ( ( - 1 , 3 ) )
self . param_vals [ self . rscale_mask ] = 0


signature:
def _poly ( self , z ) :
    """Right now legval ( z )"""

implementation:
shp = self . state . oshape . shape
zmax = float ( shp [ 0 ] )
zmin = 0.0
zmid = zmax * 0.5
coeffs = self . param_vals [ self . rscale_mask ] . copy ( )
if coeffs . size == 0 :
    ans = 0 * z
else :
    ans = np . polynomial . legendre . legval ( ( z - zmid ) / zmid , self . param_vals [ self . rscale_mask ] )
return ans


signature:
def update ( self , param_vals ) :
    """Updates all the parameters of the state + rscale ( z )"""

implementation:
self . update_rscl_x_params ( param_vals [ self . rscale_mask ] )
self . state . update ( self . param_names , param_vals [ self . globals_mask ] )
self . param_vals [ : ] = param_vals . copy ( )
if np . any ( np . isnan ( self . state . residuals ) ) :
    raise FloatingPointError ( 'state update caused nans in residuals' )


signature:
def reset ( self , ** kwargs ) :
    """Resets the aug _ state and the LMEngine"""

implementation:
self . aug_state . reset ( )
super ( LMAugmentedState , self ) . reset ( ** kwargs )


signature:
def get_shares ( self ) :
    '''Returns an object with a the numbers of shares a link has had using
 Buffer .
 www will be stripped , but other subdomains will not .'''

implementation:
self . shares = self . api . get ( url = PATHS [ 'GET_SHARES' ] % self . url ) [ 'shares' ]
return self . shares


signature:
def sample ( field , inds = None , slicer = None , flat = True ) :
    """Take a sample from a field given flat indices or a shaped slice
 Parameters
 inds : list of indices
 One dimensional ( raveled ) indices to return from the field
 slicer : slice object
 A shaped ( 3D ) slicer that returns a section of image
 flat : boolean
 Whether to flatten the sampled item before returning"""

implementation:
if inds is not None :
    out = field . ravel ( ) [ inds ]
elif slicer is not None :
    out = field [ slicer ] . ravel ( )
else :
    out = field
if flat :
    return out . ravel ( )
return out


signature:
def save ( state , filename = None , desc = '' , extra = None ) :
    """Save the current state with extra information ( for example samples and LL
 from the optimization procedure ) .
 Parameters
 state : peri . states . ImageState
 the state object which to save
 filename : string
 if provided , will override the default that is constructed based on
 the state ' s raw image file . If there is no filename and the state has
 a RawImage , the it is saved to RawImage . filename + " - peri - save . pkl "
 desc : string
 if provided , will augment the default filename to be
 RawImage . filename + ' - peri - ' + desc + ' . pkl '
 extra : list of pickleable objects
 if provided , will be saved with the state"""

implementation:
if isinstance ( state . image , util . RawImage ) :
    desc = desc or 'save'
    filename = filename or state . image . filename + '-peri-' + desc + '.pkl'
else :
    if not filename :
        raise AttributeError ( "Must provide filename since RawImage is not used" )
if extra is None :
    save = state
else :
    save = [ state ] + extra
if os . path . exists ( filename ) :
    ff = "{}-tmp-for-copy" . format ( filename )
    if os . path . exists ( ff ) :
        os . remove ( ff )
    os . rename ( filename , ff )
pickle . dump ( save , open ( filename , 'wb' ) , protocol = 2 )


signature:
def load ( filename ) :
    """Load the state from the given file , moving to the file ' s directory during
 load ( temporarily , moving back after loaded )
 Parameters
 filename : string
 name of the file to open , should be a . pkl file"""

implementation:
path , name = os . path . split ( filename )
path = path or '.'
with util . indir ( path ) :
    return pickle . load ( open ( name , 'rb' ) )


signature:
def error ( self ) :
    """Class property : Sum of the squared errors ,
 : math : ` E = \ sum _ i ( D _ i - M _ i ( \\ theta ) ) ^ 2 `"""

implementation:
r = self . residuals . ravel ( )
return np . dot ( r , r )


signature:
def loglikelihood ( self ) :
    """Class property : loglikelihood calculated by the model error ,
 : math : ` \\ mathcal { L } = - \\ frac { 1 } { 2 } \\ sum \\ left [
 \\ left ( \\ frac { D _ i - M _ i ( \\ theta ) } { \ sigma } \\ right ) ^ 2
 + \\ log { ( 2 \ pi \ sigma ^ 2 ) } \\ right ] `"""

implementation:
sig = self . hyper_parameters . get_values ( 'sigma' )
err = self . error
N = np . size ( self . data )
return - 0.5 * err / sig ** 2 - np . log ( np . sqrt ( 2 * np . pi ) * sig ) * N


signature:
def update ( self , params , values ) :
    """Update a single parameter or group of parameters ` ` params ` `
 with ` ` values ` ` .
 Parameters
 params : string or list of strings
 Parameter names which to update
 value : number or list of numbers
 Values of those parameters which to update"""

implementation: return super ( State , self ) . update ( params , values )


signature:
def push_update ( self , params , values ) :
    """Perform a parameter update and keep track of the change on the state .
 Same call structure as : func : ` peri . states . States . update `"""

implementation:
curr = self . get_values ( params )
self . stack . append ( ( params , curr ) )
self . update ( params , values )


signature:
def pop_update ( self ) :
    """Pop the last update from the stack push by
 : func : ` peri . states . States . push _ update ` by undoing the chnage last
 performed ."""

implementation:
params , values = self . stack . pop ( )
self . update ( params , values )


signature:
def temp_update ( self , params , values ) :
    """Context manager to temporarily perform a parameter update ( by using the
 stack structure ) . To use :
 with state . temp _ update ( params , values ) :
 # measure the cost or something
 state . error"""

implementation:
self . push_update ( params , values )
yield
self . pop_update ( )


signature:
def _grad_one_param ( self , funct , p , dl = 2e-5 , rts = False , nout = 1 , ** kwargs ) :
    """Gradient of ` func ` wrt a single parameter ` p ` . ( see _ graddoc )"""

implementation:
vals = self . get_values ( p )
f0 = funct ( ** kwargs )
self . update ( p , vals + dl )
f1 = funct ( ** kwargs )
if rts :
    self . update ( p , vals )
if nout == 1 :
    return ( f1 - f0 ) / dl
else :
    return [ ( f1 [ i ] - f0 [ i ] ) / dl for i in range ( nout ) ]


signature:
def _hess_two_param ( self , funct , p0 , p1 , dl = 2e-5 , rts = False , ** kwargs ) :
    """Hessian of ` func ` wrt two parameters ` p0 ` and ` p1 ` . ( see _ graddoc )"""

implementation:
vals0 = self . get_values ( p0 )
vals1 = self . get_values ( p1 )
f00 = funct ( ** kwargs )
self . update ( p0 , vals0 + dl )
f10 = funct ( ** kwargs )
self . update ( p1 , vals1 + dl )
f11 = funct ( ** kwargs )
self . update ( p0 , vals0 )
f01 = funct ( ** kwargs )
if rts :
    self . update ( p0 , vals0 )
    self . update ( p1 , vals1 )
return ( f11 - f10 - f01 + f00 ) / ( dl ** 2 )


signature:
def _grad ( self , funct , params = None , dl = 2e-5 , rts = False , nout = 1 , out = None , ** kwargs ) :
    """Gradient of ` func ` wrt a set of parameters params . ( see _ graddoc )"""

implementation:
if params is None :
    params = self . param_all ( )
ps = util . listify ( params )
f0 = funct ( ** kwargs )
# get the shape of the entire gradient to return and make an array
calc_shape = ( lambda ar : ( len ( ps ) , ) + ( ar . shape if isinstance ( ar , np . ndarray ) else ( 1 , ) ) )
if out is not None :
    grad = out
    # reference
elif nout == 1 :
    shape = calc_shape ( f0 )
    grad = np . zeros ( shape )
    # must be preallocated for mem reasons
else :
    shape = [ calc_shape ( f0 [ i ] ) for i in range ( nout ) ]
    grad = [ np . zeros ( shp ) for shp in shape ]
for i , p in enumerate ( ps ) :
    if nout == 1 :
        grad [ i ] = self . _grad_one_param ( funct , p , dl = dl , rts = rts , nout = nout , ** kwargs )
    else :
        stuff = self . _grad_one_param ( funct , p , dl = dl , rts = rts , nout = nout , ** kwargs )
        for a in range ( nout ) :
            grad [ a ] [ i ] = stuff [ a ]
return grad


signature:
def _jtj ( self , funct , params = None , dl = 2e-5 , rts = False , ** kwargs ) :
    """jTj of a ` func ` wrt to parmaeters ` params ` . ( see _ graddoc )"""

implementation:
grad = self . _grad ( funct = funct , params = params , dl = dl , rts = rts , ** kwargs )
return np . dot ( grad , grad . T )


signature:
def _hess ( self , funct , params = None , dl = 2e-5 , rts = False , ** kwargs ) :
    """Hessian of a ` func ` wrt to parmaeters ` params ` . ( see _ graddoc )"""

implementation:
if params is None :
    params = self . param_all ( )
ps = util . listify ( params )
f0 = funct ( ** kwargs )
# get the shape of the entire hessian , allocate an array
shape = f0 . shape if isinstance ( f0 , np . ndarray ) else ( 1 , )
shape = ( len ( ps ) , len ( ps ) ) + shape
hess = np . zeros ( shape )
for i , pi in enumerate ( ps ) :
    for j , pj in enumerate ( ps [ i : ] ) :
        J = j + i
        thess = self . _hess_two_param ( funct , pi , pj , dl = dl , rts = rts , ** kwargs )
        hess [ i ] [ J ] = thess
        hess [ J ] [ i ] = thess
return np . squeeze ( hess )


signature:
def crb ( self , params = None , * args , ** kwargs ) :
    """Calculate the diagonal elements of the minimum covariance of the model
 with respect to parameters params . ` ` * args ` ` and ` ` * * kwargs ` ` go to
 ` ` fisherinformation ` ` ."""

implementation:
fish = self . fisherinformation ( params = params , * args , ** kwargs )
return np . sqrt ( np . diag ( np . linalg . inv ( fish ) ) ) * self . sigma


signature:
def set_model ( self , mdl ) :
    """Setup the image model formation equation and corresponding objects into
 their various objects . ` mdl ` is a ` peri . models . Model ` object"""

implementation:
self . mdl = mdl
self . mdl . check_inputs ( self . comps )
for c in self . comps :
    setattr ( self , '_comp_' + c . category , c )


signature:
def set_image ( self , image ) :
    """Update the current comparison ( real ) image"""

implementation:
if isinstance ( image , np . ndarray ) :
    image = util . Image ( image )
if isinstance ( image , util . NullImage ) :
    self . model_as_data = True
else :
    self . model_as_data = False
self . image = image
self . _data = self . image . get_padded_image ( self . pad )
# set up various slicers and Tiles associated with the image and pad
self . oshape = util . Tile ( self . _data . shape )
self . ishape = self . oshape . pad ( - self . pad )
self . inner = self . ishape . slicer
for c in self . comps :
    c . set_shape ( self . oshape , self . ishape )
self . _model = np . zeros ( self . _data . shape , dtype = np . float64 )
self . _residuals = np . zeros ( self . _data . shape , dtype = np . float64 )
self . calculate_model ( )


signature:
def model_to_data ( self , sigma = 0.0 ) :
    """Switch out the data for the model ' s recreation of the data ."""

implementation:
im = self . model . copy ( )
im += sigma * np . random . randn ( * im . shape )
self . set_image ( util . NullImage ( image = im ) )


signature:
def get_update_io_tiles ( self , params , values ) :
    """Get the tiles corresponding to a particular section of image needed to
 be updated . Inputs are the parameters and values . Returned is the
 padded tile , inner tile , and slicer to go between , but accounting for
 wrap with the edge of the image as necessary ."""

implementation:
# get the affected area of the model image
otile = self . get_update_tile ( params , values )
if otile is None :
    return [ None ] * 3
ptile = self . get_padding_size ( otile ) or util . Tile ( 0 , dim = otile . dim )
otile = util . Tile . intersection ( otile , self . oshape )
if ( otile . shape <= 0 ) . any ( ) :
    raise UpdateError ( "update triggered invalid tile size" )
if ( ptile . shape < 0 ) . any ( ) or ( ptile . shape > self . oshape . shape ) . any ( ) :
    raise UpdateError ( "update triggered invalid padding tile size" )
# now remove the part of the tile that is outside the image and pad the
# interior part with that overhang . reflect the necessary padding back
# into the image itself for the outer slice which we will call outer
outer = otile . pad ( ( ptile . shape + 1 ) // 2 )
inner , outer = outer . reflect_overhang ( self . oshape )
iotile = inner . translate ( - outer . l )
outer = util . Tile . intersection ( outer , self . oshape )
inner = util . Tile . intersection ( inner , self . oshape )
return outer , inner , iotile


signature:
def update ( self , params , values ) :
    """Actually perform an image ( etc ) update based on a set of params and
 values . These parameter can be any present in the components in any
 number . If there is only one component affected then difference image
 updates will be employed ."""

implementation:
# FIXME needs to update priors
comps = self . affected_components ( params )
if len ( comps ) == 0 :
    return False
# get the affected area of the model image
otile , itile , iotile = self . get_update_io_tiles ( params , values )
if otile is None :
    return False
# have all components update their tiles
self . set_tile ( otile )
oldmodel = self . _model [ itile . slicer ] . copy ( )
# here we diverge depending if there is only one component update
# ( so that we may calculate a variation / difference image ) or if many
# parameters are being update ( should just update the whole model ) .
if len ( comps ) == 1 and self . mdl . get_difference_model ( comps [ 0 ] . category ) :
    comp = comps [ 0 ]
    model0 = copy . deepcopy ( comp . get ( ) )
    super ( ImageState , self ) . update ( params , values )
    model1 = copy . deepcopy ( comp . get ( ) )
    diff = model1 - model0
    diff = self . mdl . evaluate ( self . comps , 'get' , diffmap = { comp . category : diff } )
    if isinstance ( model0 , ( float , int ) ) :
        self . _model [ itile . slicer ] += diff
    else :
        self . _model [ itile . slicer ] += diff [ iotile . slicer ]
else :
    super ( ImageState , self ) . update ( params , values )
    # allow the model to be evaluated using our components
    diff = self . mdl . evaluate ( self . comps , 'get' )
    self . _model [ itile . slicer ] = diff [ iotile . slicer ]
newmodel = self . _model [ itile . slicer ] . copy ( )
# use the model image update to modify other class variables which
# are hard to compute globally for small local updates
self . update_from_model_change ( oldmodel , newmodel , itile )
return True


signature:
def get ( self , name ) :
    """Return component by category name"""

implementation:
for c in self . comps :
    if c . category == name :
        return c
return None


signature:
def _calc_loglikelihood ( self , model = None , tile = None ) :
    """Allows for fast local updates of log - likelihood"""

implementation:
if model is None :
    res = self . residuals
else :
    res = model - self . _data [ tile . slicer ]
sig , isig = self . sigma , 1.0 / self . sigma
nlogs = - np . log ( np . sqrt ( 2 * np . pi ) * sig ) * res . size
return - 0.5 * isig * isig * np . dot ( res . flat , res . flat ) + nlogs


signature:
def update_from_model_change ( self , oldmodel , newmodel , tile ) :
    """Update various internal variables from a model update from oldmodel to
 newmodel for the tile ` tile `"""

implementation:
self . _loglikelihood -= self . _calc_loglikelihood ( oldmodel , tile = tile )
self . _loglikelihood += self . _calc_loglikelihood ( newmodel , tile = tile )
self . _residuals [ tile . slicer ] = self . _data [ tile . slicer ] - newmodel


signature:
def set_mem_level ( self , mem_level = 'hi' ) :
    """Sets the memory usage level of the state .
 Parameters
 mem _ level : string
 Can be set to one of :
 * hi : all mem ' s are np . float64
 * med - hi : image , platonic are float32 , rest are float64
 * med : all mem ' s are float32
 * med - lo : image , platonic are float16 , rest float32
 * lo : all are float16 , which is bad for accuracy .
 Notes
 Right now the PSF is not affected by the mem - level changes , which is
 OK for mem but it means that self . _ model , self . _ residuals are always
 float64 , which can be a chunk of mem ."""

implementation:
# A little thing to parse strings for convenience :
key = '' . join ( [ c if c in 'mlh' else '' for c in mem_level ] )
if key not in [ 'h' , 'mh' , 'm' , 'ml' , 'm' , 'l' ] :
    raise ValueError ( 'mem_level must be one of hi, med-hi, med, med-lo, lo.' )
mem_levels = { 'h' : [ np . float64 , np . float64 ] , 'mh' : [ np . float64 , np . float32 ] , 'm' : [ np . float32 , np . float32 ] , 'ml' : [ np . float32 , np . float16 ] , 'l' : [ np . float16 , np . float16 ] }
hi_lvl , lo_lvl = mem_levels [ key ]
cat_lvls = { 'obj' : lo_lvl , 'ilm' : hi_lvl , 'bkg' : lo_lvl }
# no psf . . .
self . image . float_precision = hi_lvl
self . image . image = self . image . image . astype ( lo_lvl )
self . set_image ( self . image )
for cat in cat_lvls . keys ( ) :
    obj = self . get ( cat )
    # check if it ' s a component collection
    if hasattr ( obj , 'comps' ) :
        for c in obj . comps :
            c . float_precision = lo_lvl
    else :
        obj . float_precision = lo_lvl
self . _model = self . _model . astype ( hi_lvl )
self . _residuals = self . _model . astype ( hi_lvl )
self . reset ( )


signature:
def scramble_positions ( p , delete_frac = 0.1 ) :
    """randomly deletes particles and adds 1 - px noise for a realistic
 initial featuring guess"""

implementation:
probs = [ 1 - delete_frac , delete_frac ]
m = np . random . choice ( [ True , False ] , p . shape [ 0 ] , p = probs )
jumble = np . random . randn ( m . sum ( ) , 3 )
return p [ m ] + jumble


signature:
def create_img ( ) :
    """Creates an image , as a ` peri . util . Image ` , which is similar
 to the image in the tutorial"""

implementation:
# 1 . particles + coverslip
rad = 0.5 * np . random . randn ( POS . shape [ 0 ] ) + 4.5
# 4.5 + - 0.5 px particles
part = objs . PlatonicSpheresCollection ( POS , rad , zscale = 0.89 )
slab = objs . Slab ( zpos = 4.92 , angles = ( - 4.7e-3 , - 7.3e-4 ) )
objects = comp . ComponentCollection ( [ part , slab ] , category = 'obj' )
# 2 . psf , ilm
p = exactpsf . FixedSSChebLinePSF ( kfki = 1.07 , zslab = - 29.3 , alpha = 1.17 , n2n1 = 0.98 , sigkf = - 0.33 , zscale = 0.89 , laser_wavelength = 0.45 )
i = ilms . BarnesStreakLegPoly2P1D ( npts = ( 16 , 10 , 8 , 4 ) , zorder = 8 )
b = ilms . LegendrePoly2P1D ( order = ( 7 , 2 , 2 ) , category = 'bkg' )
off = comp . GlobalScalar ( name = 'offset' , value = - 2.11 )
mdl = models . ConfocalImageModel ( )
st = states . ImageState ( util . NullImage ( shape = [ 48 , 64 , 64 ] ) , [ objects , p , i , b , off ] , mdl = mdl , model_as_data = True )
b . update ( b . params , BKGVALS )
i . update ( i . params , ILMVALS )
im = st . model + np . random . randn ( * st . model . shape ) * 0.03
return util . Image ( im )


signature:
def get_values ( self , params ) :
    """Get the value of a list or single parameter .
 Parameters
 params : string , list of string
 name of parameters which to retrieve"""

implementation: return util . delistify ( [ self . param_dict [ p ] for p in util . listify ( params ) ] , params )


signature:
def set_values ( self , params , values ) :
    """Directly set the values corresponding to certain parameters .
 This does not necessarily trigger and update of the calculation ,
 See also
 : func : ` ~ peri . comp . comp . ParameterGroup . update ` : full update func"""

implementation:
for p , v in zip ( util . listify ( params ) , util . listify ( values ) ) :
    self . param_dict [ p ] = v


signature:
def set_shape ( self , shape , inner ) :
    """Set the overall shape of the calculation area . The total shape of that
 the calculation can possibly occupy , in pixels . The second , inner , is
 the region of interest within the image ."""

implementation:
if self . shape != shape or self . inner != inner :
    self . shape = shape
    self . inner = inner
    self . initialize ( )


signature:
def trigger_update ( self , params , values ) :
    """Notify parent of a parameter change"""

implementation:
if self . _parent :
    self . _parent . trigger_update ( params , values )
else :
    self . update ( params , values )


signature:
def split_params ( self , params , values = None ) :
    """Split params , values into groups that correspond to the ordering in
 self . comps . For example , given a sphere collection and slab : :
 ( spheres ) [ pos rad etc ] [ pos val , rad val , etc ]
 ( slab ) [ slab params ] [ slab vals ]"""

implementation:
pc , vc = [ ] , [ ]
returnvalues = values is not None
if values is None :
    values = [ 0 ] * len ( util . listify ( params ) )
for c in self . comps :
    tp , tv = [ ] , [ ]
    for p , v in zip ( util . listify ( params ) , util . listify ( values ) ) :
        if not p in self . lmap :
            raise NotAParameterError ( "%r does not belong to %r" % ( p , self ) )
        if c in self . pmap [ p ] :
            tp . append ( p )
            tv . append ( v )
    pc . append ( tp )
    vc . append ( tv )
if returnvalues :
    return pc , vc
return pc


signature:
def get ( self ) :
    """Combine the fields from all components"""

implementation:
fields = [ c . get ( ) for c in self . comps ]
return self . field_reduce_func ( fields )


signature:
def set_shape ( self , shape , inner ) :
    """Set the shape for all components"""

implementation:
for c in self . comps :
    c . set_shape ( shape , inner )


signature:
def sync_params ( self ) :
    """Ensure that shared parameters are the same value everywhere"""

implementation:
def _normalize ( comps , param ) :
    vals = [ c . get_values ( param ) for c in comps ]
    diff = any ( [ vals [ i ] != vals [ i + 1 ] for i in range ( len ( vals ) - 1 ) ] )
    if diff :
        for c in comps :
            c . set_values ( param , vals [ 0 ] )
for param , comps in iteritems ( self . lmap ) :
    if isinstance ( comps , list ) and len ( comps ) > 1 :
        _normalize ( comps , param )


signature:
def setup_passthroughs ( self ) :
    """Inherit some functions from the components that we own . In particular ,
 let ' s grab all functions that begin with ` param _ ` so the super class
 knows how to get parameter groups . Also , take anything that is listed
 under Component . exports and rename with the category type , i . e . ,
 SphereCollection . add _ particle - > Component . obj _ add _ particle"""

implementation:
self . _nopickle = [ ]
for c in self . comps : # take all member functions that start with ' param _ '
    funcs = inspect . getmembers ( c , predicate = inspect . ismethod )
    for func in funcs :
        if func [ 0 ] . startswith ( 'param_' ) :
            setattr ( self , func [ 0 ] , func [ 1 ] )
            self . _nopickle . append ( func [ 0 ] )
    # add everything from exports
    funcs = c . exports ( )
    for func in funcs :
        newname = c . category + '_' + func . __func__ . __name__
        setattr ( self , newname , func )
        self . _nopickle . append ( newname )


signature:
def get_conf_filename ( ) :
    """The configuration file either lives in ~ / . peri . json or is specified on the
 command line via the environment variables PERI _ CONF _ FILE"""

implementation:
default = os . path . join ( os . path . expanduser ( "~" ) , ".peri.json" )
return os . environ . get ( 'PERI_CONF_FILE' , default )


signature:
def read_environment ( ) :
    """Read all environment variables to see if they contain PERI"""

implementation:
out = { }
for k , v in iteritems ( os . environ ) :
    if transform ( k ) in default_conf :
        out [ transform ( k ) ] = v
return out


signature:
def load_conf ( ) :
    """Load the configuration with the priority :
 1 . environment variables
 2 . configuration file
 3 . defaults here ( default _ conf )"""

implementation:
try :
    conf = copy . copy ( default_conf )
    conf . update ( json . load ( open ( get_conf_filename ( ) ) ) )
    conf . update ( read_environment ( ) )
    return conf
except IOError as e :
    create_default_conf ( )
    return load_conf ( )


signature:
def get_group_name ( id_group ) :
    """Used for breadcrumb dynamic _ list _ constructor ."""

implementation:
group = Group . query . get ( id_group )
if group is not None :
    return group . name


signature:
def index ( ) :
    """List all user memberships ."""

implementation:
page = request . args . get ( 'page' , 1 , type = int )
per_page = request . args . get ( 'per_page' , 5 , type = int )
q = request . args . get ( 'q' , '' )
groups = Group . query_by_user ( current_user , eager = True )
if q :
    groups = Group . search ( groups , q )
groups = groups . paginate ( page , per_page = per_page )
requests = Membership . query_requests ( current_user ) . count ( )
invitations = Membership . query_invitations ( current_user ) . count ( )
return render_template ( 'invenio_groups/index.html' , groups = groups , requests = requests , invitations = invitations , page = page , per_page = per_page , q = q )


signature:
def requests ( ) :
    """List all pending memberships , listed only for group admins ."""

implementation:
page = request . args . get ( 'page' , 1 , type = int )
per_page = request . args . get ( 'per_page' , 5 , type = int )
memberships = Membership . query_requests ( current_user , eager = True ) . all ( )
return render_template ( 'invenio_groups/pending.html' , memberships = memberships , requests = True , page = page , per_page = per_page , )


signature:
def invitations ( ) :
    """List all user pending memberships ."""

implementation:
page = request . args . get ( 'page' , 1 , type = int )
per_page = request . args . get ( 'per_page' , 5 , type = int )
memberships = Membership . query_invitations ( current_user , eager = True ) . all ( )
return render_template ( 'invenio_groups/pending.html' , memberships = memberships , page = page , per_page = per_page , )


signature:
def new ( ) :
    """Create new group ."""

implementation:
form = GroupForm ( request . form )
if form . validate_on_submit ( ) :
    try :
        group = Group . create ( admins = [ current_user ] , ** form . data )
        flash ( _ ( 'Group "%(name)s" created' , name = group . name ) , 'success' )
        return redirect ( url_for ( ".index" ) )
    except IntegrityError :
        flash ( _ ( 'Group creation failure' ) , 'error' )
return render_template ( "invenio_groups/new.html" , form = form , )


signature:
def manage ( group_id ) :
    """Manage your group ."""

implementation:
group = Group . query . get_or_404 ( group_id )
form = GroupForm ( request . form , obj = group )
if form . validate_on_submit ( ) :
    if group . can_edit ( current_user ) :
        try :
            group . update ( ** form . data )
            flash ( _ ( 'Group "%(name)s" was updated' , name = group . name ) , 'success' )
        except Exception as e :
            flash ( str ( e ) , 'error' )
            return render_template ( "invenio_groups/new.html" , form = form , group = group , )
    else :
        flash ( _ ( 'You cannot edit group %(group_name)s' , group_name = group . name ) , 'error' )
return render_template ( "invenio_groups/new.html" , form = form , group = group , )


signature:
def delete ( group_id ) :
    """Delete group ."""

implementation:
group = Group . query . get_or_404 ( group_id )
if group . can_edit ( current_user ) :
    try :
        group . delete ( )
    except Exception as e :
        flash ( str ( e ) , "error" )
        return redirect ( url_for ( ".index" ) )
    flash ( _ ( 'Successfully removed group "%(group_name)s"' , group_name = group . name ) , 'success' )
    return redirect ( url_for ( ".index" ) )
flash ( _ ( 'You cannot delete the group %(group_name)s' , group_name = group . name ) , 'error' )
return redirect ( url_for ( ".index" ) )


signature:
def members ( group_id ) :
    """List user group members ."""

implementation:
page = request . args . get ( 'page' , 1 , type = int )
per_page = request . args . get ( 'per_page' , 5 , type = int )
q = request . args . get ( 'q' , '' )
s = request . args . get ( 's' , '' )
group = Group . query . get_or_404 ( group_id )
if group . can_see_members ( current_user ) :
    members = Membership . query_by_group ( group_id , with_invitations = True )
    if q :
        members = Membership . search ( members , q )
    if s :
        members = Membership . order ( members , Membership . state , s )
    members = members . paginate ( page , per_page = per_page )
    return render_template ( "invenio_groups/members.html" , group = group , members = members , page = page , per_page = per_page , q = q , s = s , )
flash ( _ ( 'You are not allowed to see members of this group %(group_name)s.' , group_name = group . name ) , 'error' )
return redirect ( url_for ( '.index' ) )


signature:
def leave ( group_id ) :
    """Leave group ."""

implementation:
group = Group . query . get_or_404 ( group_id )
if group . can_leave ( current_user ) :
    try :
        group . remove_member ( current_user )
    except Exception as e :
        flash ( str ( e ) , "error" )
        return redirect ( url_for ( '.index' ) )
    flash ( _ ( 'You have successfully left %(group_name)s group.' , group_name = group . name ) , 'success' )
    return redirect ( url_for ( '.index' ) )
flash ( _ ( 'You cannot leave the group %(group_name)s' , group_name = group . name ) , 'error' )
return redirect ( url_for ( '.index' ) )


signature:
def approve ( group_id , user_id ) :
    """Approve a user ."""

implementation:
membership = Membership . query . get_or_404 ( ( user_id , group_id ) )
group = membership . group
if group . can_edit ( current_user ) :
    try :
        membership . accept ( )
    except Exception as e :
        flash ( str ( e ) , 'error' )
        return redirect ( url_for ( '.requests' , group_id = membership . group . id ) )
    flash ( _ ( '%(user)s accepted to %(name)s group.' , user = membership . user . email , name = membership . group . name ) , 'success' )
    return redirect ( url_for ( '.requests' , group_id = membership . group . id ) )
flash ( _ ( 'You cannot approve memberships for the group %(group_name)s' , group_name = group . name ) , 'error' )
return redirect ( url_for ( '.index' ) )


signature:
def remove ( group_id , user_id ) :
    """Remove user from a group ."""

implementation:
group = Group . query . get_or_404 ( group_id )
user = User . query . get_or_404 ( user_id )
if group . can_edit ( current_user ) :
    try :
        group . remove_member ( user )
    except Exception as e :
        flash ( str ( e ) , "error" )
        return redirect ( urlparse ( request . referrer ) . path )
    flash ( _ ( 'User %(user_email)s was removed from %(group_name)s group.' , user_email = user . email , group_name = group . name ) , 'success' )
    return redirect ( urlparse ( request . referrer ) . path )
flash ( _ ( 'You cannot delete users of the group %(group_name)s' , group_name = group . name ) , 'error' )
return redirect ( url_for ( '.index' ) )


signature:
def accept ( group_id ) :
    """Accpet pending invitation ."""

implementation:
membership = Membership . query . get_or_404 ( ( current_user . get_id ( ) , group_id ) )
# no permission check , because they are checked during Memberships creating
try :
    membership . accept ( )
except Exception as e :
    flash ( str ( e ) , 'error' )
    return redirect ( url_for ( '.invitations' , group_id = membership . group . id ) )
flash ( _ ( 'You are now part of %(name)s group.' , user = membership . user . email , name = membership . group . name ) , 'success' )
return redirect ( url_for ( '.invitations' , group_id = membership . group . id ) )


signature:
def new_member ( group_id ) :
    """Add ( invite ) new member ."""

implementation:
group = Group . query . get_or_404 ( group_id )
if group . can_invite_others ( current_user ) :
    form = NewMemberForm ( )
    if form . validate_on_submit ( ) :
        emails = filter ( None , form . data [ 'emails' ] . splitlines ( ) )
        group . invite_by_emails ( emails )
        flash ( _ ( 'Requests sent!' ) , 'success' )
        return redirect ( url_for ( '.members' , group_id = group . id ) )
    return render_template ( "invenio_groups/new_member.html" , group = group , form = form )
flash ( _ ( 'You cannot invite users or yourself (i.e. join) to the group ' '%(group_name)s' , group_name = group . name ) , 'error' )
return redirect ( url_for ( '.index' ) )