([ang-step],[ang-limit],[radius],[alignsh],[prj-radius],[iter],[grp],[toobig],[maxspfreq]) ; Main refinement group loop
;
; SOURCE: grploop.pam    Original                            ArDean Leith  Nov 2000
;                        %degredations                       ArDean Leith  Feb 2005
;                        []                                  ArDean Leith  Dec 2005
;                        More stacks & 'RT SQ' selection     ArDean Leith  Dec 2006
;                        'AP SH' use                         ArDean Leith  Feb 2007
;                        tmp dir purge                       ArDean Leith  Mar 2008
;
; MASTER COPY: /net/bali/usr1/spider/docs/techs/recon/newprogs/
; 
; PURPOSE: Main refinement group loop.  Runs for each defocus group on each iteration
;
; CALLED FROM: pub_refine_start or   
;              refine.pam  
;
; INPUT REGISTERS:
;    [ang-step]             Angular steps           (varies with iter)
;    [ang-limit]            Restrict angular search (varies with iter)
;    [radius]               Radius of the structure (pixels)
;    [alignsh]              Shift allowed is +-[alignsh]
;    [prj-radius]           Radius of the structure used in projection
;    [iter]                 Alignment step iteration counter  (varies with iter)
;    [grp]                  Defocus group                     (varies with group)
;
; OUTPUT REGISTERS:
;    [toobig]               % of images whose proj. angle moved by > 1.5*[ang-step]
;
;  '##' denotes iteration,  '##+' denotes next iteration, and '***' denotes group
; INPUT FILES (SET IN: refine_settings  refine_settings.pam.pam):
;    [order_select]         input/order_select           Image ID doc. fil
;    [sorted_order_select]  input/order_select_sort      Image ID doc. fil
;    [sel_particles]        input/select_***             Particle selection file
;    [iter_vft]             final/vft##                  Current filtered volume 
;    [temp_ctf_file]        work/ctf***                  CTF corrected volume from prepare
;    [aligned_images]       input/dala##_***             Aligned stacked image files
;    [group_align]          final/align##+_***           Alignment parameter doc file     
;    [unaligned_images]     input/data***                Unaligned stacked image files  
;    [iter_ang_voea]        final/angvoea_##             Ref. angles doc. file        
;    [iter_select_voea]     final/selvoea_##             Ref. angles selection doc. file    
;
; OUTPUT FILES (SET IN: refine_settings refine_settings.pam.pam):
;    [next_group_align]     final/align##_***            Alignment parameter doc file
;    [next_aligned_images]  work/dala##+_***@            Aligned stacked image files
;    [next_group_vol]       work/defgrp***/vol##         Reconstructed volume
;    [next_group_dres]      final/defgrp***/dres##+      Resolution curve doc. file
;    [next_group_vol]_odd   work/defgrp***/vol_##+_odd   (created & deleted)
;    [next_group_vol]_even  work/defgrp***/vol_##+_even  (created & deleted)
;    [temp_ref_projs]       /tmp/refproj##_***           (created & deleted)
; 
; PROCEDURES CALLED: saveresp
;
; INLINE BUFFERS USED: _1, _4, _8

 VM
 echo " In grploop, Iteration: {**[iter]}  Group: {***[grp]}"
 VM
 date '+ Time: %x  %X'

 [radius1]=5.0                             ; First radius for AP (Can alter this)
 [ang-change-thresh]=1.5*[ang-step]        ; Convergence criterion ang. distance limit.
 [next-iter]=[iter]+1
  
 UD N [num-refs]                           ; Get number of reference images used
 [iter_select_voea]                        ; Ref. angles selection doc. file  (input)

 UD N,[numgrps]                            ; Find number of defocus groups
 [sorted_order_select]                     ;                                  (input)

 ; Get number of images in this group
 DO LB1 [i]=1,[numgrps]                    ; Loop over defocus groups  -----------
    ;         GROUP, PART.               
    UD IC [i], [grpt], [numparts]          ; Get number of particles 
    [sorted_order_select]                  ; Doc file listing groups (input)    
    IF ([grpt].EQ.[grp])GOTO LB2
 LB1
 LB2                                       ; Got number of particles in this group

 UD ICE
 [sorted_order_select]

 ; Multiply Fourier of current volume by CTF file for this group
 MU                                        ; Multiply
 [iter_vft]                                ; Fourier of current volume         (input}
 [temp_ctf_file]                           ; CTF  file     (work/ctf{grp})     (input)
 _1                                        ; _1 created here                   (output)
 *                                         ; End multiplications
 
 FT                                        ; Fourier back transform
 _1                                        ; Fourier of CTF corrected volume   (input)
 _4                                        ; CTF corrected current volume      (output) 

 DE                                        ; _1 can be deleted here 
 _1                                        ; Fourier of CTF corrected current volume

 MY FL                                     ; Flush results file

 ; Clear tmp directory and create temp local dir., divert errors to junk
 VM
 \rm -rf /tmp/* >& /dev/null

 VM
 mkdir -p [temp_local_dir]

 VM                                        ; So anyone can delete it
 chmod 777 [temp_local_dir]

 ; Create stack holding angular reference projections from CTF corrected volume.
 PJ 3Q                                     ; Create reference projections
 _4                                        ; CTF corrected current volume     (input)
 [prj-radius]                              ; Radius of object
 [iter_select_voea]                        ; Ref. angles selection doc. file  (input) 
 [iter_ang_voea]                           ; Ref. angles doc. file            (input)
 [temp_ref_projs]@******                   ; Template for ref. projections    (output) 

 DE                                        ; _4 can be deleted here
 _4                                        ; CTF corrected current volume

 FI X [maxim]                              ; Find total number of images (not [numparts])
 [aligned_images]                          ; Input file needed                (input)
 (26)                                      ; Max. image number header position

 CP                                        ; Copy current aligned images to inline stack
 [aligned_images]                          ; Input file needed                (input)
 _8@                                       ; Current aligned images stack     (output)
 [maxim]                                   ; Number of images in stack _8@
 
 DE                                        ; Remove existing alignment doc file  
 [next_group_align]                        ; Alignment parameter doc file

 DE                                        ; Remove existing scratch file  
 work/SCR_{**[iter]}_{***[grp]}            ; Reference rings scratch file 

 MY FL                                     ; Flush results file

 [useapsh]= 0                              ; Do not use 'AP SH' (Can change this)
 IF ([iter] .GT. 1) THEN
    [useapsh] = 0                          ; Use 'AP REF'instead of 'AP SH'
 ENDIF

 ; Find reference projection matching current aligned image
 ;   (For large images change 'skip' to 2 or 3 to decrease memory)
 IF ([useapsh] .EQ. 1) THEN
    ; 'AP SH' gives more exhaustive search than 'AP REF' but is 4-5x slower
    AP SH 
    [temp_ref_projs]@******                ; Template for ref. projections    (input)
    (1-[num-refs])                         ; Ref. projections file numbers
    ([alignsh],1)                          ; Shift search range, Step size
    ([radius1],[radius],1)                 ; First, last radial ring, & skip
    [iter_ang_voea]                        ; Ref. angles file                 (input)
    _8@******                              ; Template for exp. aligned images (input)
    [sel_particles]                        ; Particle selection files         (input)
    [group_align]                          ; Alignment parameter doc. file    (input)
    [ang-limit],[ang-change-thresh]        ; Angular search restriction
    (1)                                    ; Check mirrored projections
    [next_group_align]                     ; Alignment parameter doc file     (output)

    ; Default for  unknown return value
    [toobig]= 100

    VM
    echo " Iteration: {**[iter]} Group: {*****[grp]} "

 ELSE

    AP REF 
    [temp_ref_projs]@******                ; Template for ref. projections    (input)
    (1-[num-refs])                         ; Ref. projections file numbers
    [alignsh]                              ; Shift search range
    ([radius1],[radius],1)                 ; First, last radial ring, & skip
    [iter_ang_voea]                        ; Ref. angles file                 (input)
    work/SCR_{**[iter]}_{***[grp]}         ; No scratch file if fits in-core
    _8@******                              ; Template for exp. aligned images (input)
    [sel_particles]                        ; Particle selection files         (input)
    [group_align]                          ; Alignment parameter doc. file    (input)
    [ang-limit],[ang-change-thresh]        ; Angular search restriction
    (1)                                    ; Check mirrored projections
    [next_group_align]                     ; Alignment parameter doc file     (output)

    ; Check size of change in proj. angle
    ;     %BIG-ANGDIF,       AVG-ANGDIF,  AVG-CCROT, %WORSE,  AVG-WORSE,   AVG-BETTER      
    UD -2,[percent-greater],[AVG-ANGDIF],[avg-ccrot],[degred],[avg-degred],[avg-impr]
    [next_group_align]                        ; Alignment parameter doc file      (input)
    UD E                                      ; Close doc file access

    ; Return % of image whose proj. angle moved by > 1.5*[ang-step] criterion
    [toobig]=[percent-greater]* 100

    VM
    echo " Iteration: {**[iter]} Group: {*****[grp]}  Excessive changes: {%f8.2%[percent-greater]}% "

    VM
    echo " CCROT Degredation: {***[degred]}%   Avg. degred.: {%F8.2%[avg-degred]}  Avg. improve.: {%F8.2%[avg-impr]}"
 ENDIF

 ; Apply transformations to all original particle images
 RT SQ                                     ; Rotate & shift operation
 [unaligned_images]@******                 ; Unaligned original stacked images 
 [sel_particles]                           ; Particle selection file           (input) 
 (6,0,7,8)                                 ; Reg. numbers for angle, scale, & shift
 [next_group_align]                        ; Alignment parameter doc file      (input)
 _8@******                                 ; Current aligned images            (output)
 
 ; Store aligned projections (Do now in case BP fails)
 CP                                        ; Create hard copy of aligned images
 _8@                                       ; Aligned images in inline file     (input)
 [next_aligned_images]                     ; On disk file                      (output)

 VM                                        ; Echo for isolating stack overflow
 echo " Back projecting group: {*****[grp]}"

 MY FL

 ; Calculate new, refined volume using centered projections and 
 ; the corrected angles from align doc. file.  Creates two additional
 ; volumes from odd/even images for use in resolution calculation.

 ; (If you have large images which give problems allocating memory in 'BP 32F', 
 ;     substitute operation 'BP 3F'.  Use that operation three times (with 3 
 ;     appropriate selection files for the images to be included) to create 
 ;     the three output volumes one by one.)
 BP 32F                         ; Back Projection - 3D Fourier
 _8@******                      ; Template for current aligned images  (input)
 [sel_particles]                ; Particle selection doc. file         (input)
 [next_group_align]             ; Alignment parameter doc file         (input)
 *                              ; No symmetries  
 [next_group_vol]tmp            ; Reconstructed 3D file                (output)
 [next_group_vol]_odd           ; Reconstructed 3D file                (output)
 [next_group_vol]_even          ; Reconstructed 3D file                (output)

 DE                             ; _8 no longer needed
 _8

 MY FL                          ; Flush results file

 ; Find center of gravity of the temp. volume
 CG PH [unused],[unused],[unused],[xcg],[ycg],[zcg]  ; Center of Gravity of vol.
 [next_group_vol]tmp            ; Volume                             (input)

 SH F                           ; Shift temp volume to its center of gravity
 [next_group_vol]tmp            ; Temp. volume  file                 (input)
 [next_group_vol]               ; Final group volume                 (output)
 -[xcg],-[ycg],-[zcg]           ; Shift

 DE                             ; Delete temp volume file 
 [next_group_vol]tmp

 RF 3 [unused],[spfreq]         ; Find phase residual &  shell correl.
 [next_group_vol]_odd           ; First  volume file                (input)
 [next_group_vol]_even          ; Second volume file                (input)
 (0.5)                          ; Ring width
 (0.2,2.0)                      ; Lower and upper scale factor
 C                              ; Missing cone
 (90.0)                         ; Maximum tilt angle 
 (3)                            ; Factor for noise comparison
 [next_group_dres]              ; Resolution document file          (output)

 ; Record this groups reconstruction resolution in doc file
 @saveresp([maxspfreq],[iter],[grp],[spfreq])
 [grp_resol]                    ; Resolution summary file           (output)

 VM                             ; Remove - refproj (no longer needed)
 \rm -f [temp_ref_projs]*.$DATEXT 

 IF ([iter] .GT. 1) THEN        ; Not first iteration
    DE                          ; Remove previous group vol. 
    [group_vol]                 ;   work/vol_{**[iter]}_{***[grp]}      
    DE                          ; Remove previous odd group vol.
    [group_vol]_odd             ;   work/vol_{**[iter]}_{***[grp]}_odd
    DE                          ; Remove previous even group vol.
    [group_vol]_even            ;   work/vol_{**[iter]}_{***[grp]}_even
    DE                          ; Remove previous current aligned images
    [aligned_images_prefix]     ;   work/dala{**[iter]}_{***[grp]}
 ENDIF

 MY FL
RE
;