; RUNS CORRESPONDENCE ANALYSIS OR PRINCIPLE COMPONENT ANALYSIS

; NOTES: 
;   uses Python script eigendoc.py (jsl) & Gnuplot script ploteigen.gnu (jsl & trs)

; ---------------- Parameters ----------------

x65 = 65   ; image dimension (square)
x32 = 32   ; particle radius (pixels)
x29 = .05  ; additive constant for correspondence analysis
x27 = 9    ; number of eigenfactors to calculate

; ------------------ Inputs ------------------
fr l
[particles]ali/sar******      ; particle file-pattern
fr l
[selection_doc]listparticles  ; selection doc file

; ------------------ Outputs ------------------
fr l
[ca_dir]CA                                    ; output directory
fr l
[cas_prefix][ca_dir]/cas                      ; correspondence-analysis output prefix
fr l
[ps_eigenvalue_hist][ca_dir]/eigenvalues      ; eigenvalue histogram, encapsulated PostScript
fr l
[eigen_img][ca_dir]/eigenimg***               ; eigenimage template
fr l
[reconstituted_img][ca_dir]/reconstituted***  ; reconstituted image template
fr l
[ps_factor_map][ca_dir]/factormap             ; eigenvalue factor map, PostScript format
fr l
[sdc_doc][ca_dir]/doccorrmap                  ; eigenvalue (1-9) document file (for WEB->CorrMap)

; ------------- END BATCH HEADER -------------

vm
echo "if(! -d [ca_dir]) mkdir [ca_dir]"|csh

md
set mp
0

; make circular mask
mo  
_6       ; output mask file
x65,x65  ; dimensions
C        ; _C_ircle
x32      ; radius


; CHECK FOR SMALL/NEGATIVE PIXEL INTENSITIES

vm
echo "Checking image statistics"; date

x79=9999  ; initialize minimum

; loop through particles
do lb5 x50=1,x27
    ; get particle number (x52)
    ud ic,x50,x52
    [selection_doc]

    ; get max (x88)
    fs
    [particles]x52
    fi x88
    [particles]x52
    (8)  ; header location for minimum intensity

    ; update if necessary
    if(x88.lt.x79) then
        x79=x88  ; minimum intensity
        x49=x52  ; particle#
    endif
lb5
; end particle-loop

ud ice  ; close document
[selection_doc]

vm m
echo "Minimum intensity: {%f8.4%x79} (image #{******x49})"; 
echo ;
echo "Running correspondence analysis"; date
.

; run correspondence analysis
ca s
[particles]
[selection_doc]
_6            ; mask
x27           ; number of eigenfactors
C             ; _C_orrespondence analysis (P<-PCA)
x29           ; additive factor (remove line if using PCA)
[cas_prefix]  ; output file prefix

lb1

; plot eigenvalues
vm m
echo "Preparing eigenvalue plot"; 
python eigendoc.py [cas_prefix]_EIG.$DATEXT > eigenvalues.txt ;
sed "s/LAST/{**x27}/" ploteigen.gnu |
sed "s,EIGENVALUES,[ps_eigenvalue_hist],"  > tmp2.gnu ;
gnuplot tmp2.gnu ;
rm -f eigenvalues.txt tmp2.gnu
.


; RECONSTITUTE EIGENIMAGES

vm
echo "Reconstituting eigenimages"; date

x13=x65*2  ; double image-dimension
x66=x65+1  ; image-dimension + 1

; loop through eigenfactors
do lb9 x90=1,x27
    ; generate eigenimage
    ca sre
    [cas_prefix]
    x90  ; eigenvector
    [eigen_img]x90

    ; generate positive reconstituted image
    ca sra
    [cas_prefix]
    x90    ; eigenvector
    (0.2)  ; eigenvalue
    _1     ; OUTPUT

    ; pad image to twice the height
    pd
    _1  ; INPUT
    [reconstituted_img]x90
    x65,x13  ; xdim, ydim
    B        ; set background to _B_order

    ; generate negative reconstituted image
    ca sra
    [cas_prefix]
    x90     ; eigenvector
    (-0.2)  ; eigenvalue
    _3      ; OUTPUT

    ; insert negative reconstituted image into larger image
    in
    _3                      ; INPUT: small, negative reconstituted image
    [reconstituted_img]x90  ; INPUT (overwritten)
    (1,x66)                 ; x-,y-coords
lb9


; PREPARE FACTOR MAPS

vm
echo "Generating factor maps"; date

x26=x27-1  ; number of eigenfactor-pairs to plot

; loop through eigenfactor-pairs
do lb2 x21=1,x26
    x22=x21+1  ; second factor to plot

    ; generate factor map
    ca sm
    I        ; plot _I_mage coordinates
    [cas_prefix]  ; INPUT
    (0)      ; #horizontal patches
    x21,x22  ; factors to plot
    S        ; plot _S_ymbol
    +        ; symbol to plot
    Y        ; output to PostScript?
    (2.3)    ; #standard devations to plot
    (0)      ; no axis-flip
    [ps_factor_map]{**x21}{**x22}  ; OUTPUT
             ; TEXT SIZE FOR AXIS AND DATA
             ; X AXIS OFFSET
             ; LOWER, UPPER AXIS BOUNDS
             ; AXIS LABEL UNIT AND TICKS / LABEL
             ; LABEL NO. TO EDIT
             ; Y AXIS OFFSET
             ; LOWER, UPPER AXIS BOUNDS
             ; AXIS LABEL UNIT AND TICKS / LABEL
             ; LABEL NO. TO EDIT
lb2
; end pair-loop


; write eigenvalues to document file
sd c
[cas_prefix]
1-9  ; eigenfactor range (past nine crashes)
[sdc_doc]

vm
echo "Done"; date

en d

; Modified 2007-11-02
;   2007-11-02 (trs) -- added reconstituted images
;   2006-05-11 (trs) -- bug fix if more than nine eigenfactors
;   2005-08-23 (trs) -- factor maps adapted from jsl's casm.bat