Inverse masking code

The selection (or masking) operator $\bold K$ is implemented in mask1() .  

module mask1 { # masking operator
  logical, dimension( :), pointer :: m
#% _init( m)
#% _lop( x, y)
  if( adj) 
	   where( m) x += y
  else  #
           where( m) y += x
}

The inverting of the mask operator proceeds much as we inverted the linear interpolation operator with module invint2 . The main difference is we swap the selection operator for the linear interpolation operator. (Philosophically, selection is like binning which is like nearest-neighbor interpolation.) The module mis2 , does the job.

 

module mis2 {
  use mask1 + helicon + polydiv + cgstep_mod + solver_mod
contains
# fill in missing data by minimizing power out of a given filter.
# by helix magic works in any number of dimensions
  subroutine mis1( niter, xx, aa, known, doprec) {
    logical,                intent( in)     :: doprec
    integer,                intent( in)     :: niter
    type( filter),          intent( in)     :: aa
    logical, dimension( :), intent( in)     :: known
    real,    dimension( :), intent( in out) :: xx	# fitting variables
    real,    dimension( :), allocatable     :: dd
    logical, dimension( :), pointer         :: kk
    integer                                 :: nx
    nx = size( xx)
    if( doprec) {                          #  preconditioned
       allocate( kk( nx)); kk = known
       call mask1_init( kk)
       call polydiv_init( nx, aa)
       call solver_prec( mask1_lop, cgstep, niter= niter, x= xx, dat= xx, 
                    prec= polydiv_lop, nprec= nx, eps= 0.)
       call polydiv_close()
       deallocate( kk)
    } else {                               #  regularized
       allocate( dd( nx)); dd = 0.
       call helicon_init( aa)
       call solver( helicon_lop, cgstep, niter= niter, x= xx, dat= dd, 
                    known = known, x0= xx)
       deallocate( dd)
    }
    call cgstep_close( )
  }
}

It is instructive to compare mis2 with invint2 . Both are essentially filling empty regions consistant with prior knowledge at particular locations and minimizing energy of the filtered field. Both use the helix and can be used in N-dimensional space.