Inverse masking code

The selection (or masking) operator 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.

 

#$=head1 NAME
#$
#$mis2 - Fill in missing data
#$
#$=head1 SYNOPSIS 
#$  
#$C<call mis1(niter,xx,aa,known,doprec)>  
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$  
#$=item  niter - integer   
#$  
#$       number of iterations
#$
#$=item  xx    -real
#$
#$       fitting variable
#$
#$=item  aa    -type(filter)
#$
#$       filter to apply
#$  
#$=item  known -C<logical(:)>
#$  
#$       Known data
#$
#$=item  doprec -logical
#$
#$       Whether or not to run preconditioning
#$     
#$=back   
#$     
#$=head1 DESCRIPTION
#$     
#$fill in missing data by minimizing power out of a given filter
#$ by helix magic works in any number of dimensions
#$  
#$=head1 SEE ALSO
#$
#$L<nmis2>,L<msmis>,L<mask1>,L<helicon>,L<solver>
#$  
#$=head1 LIBRARY
#$  
#$B<geef90>
#$
#$=cut
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.