; 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]../Particles/flt/flt  ; particle file-pattern, ends in six digits
fr l
[selection_doc]select/sel001     ; 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, will end in three digits
fr l
[ps_factor_map][ca_dir]/factormap         ; eigenvalue factor map, PostScript format
fr l
[sdc_doc][ca_dir]/doccorrmap              ; eigenvalue 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]{***x80}

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

; loop through eigenfactors
do lb9 x90=1,x27
    ca sre
    [cas_prefix]
    x90  ; eigenvector(s)
    [eigen_img]{***x90}
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-x27  ; eigenfactor range
[sdc_doc]

vm
echo "\nDone"; date

en d

; Modified 2005-09-02
;   2005-08-23 (trs) -- factor maps adapted from jsl's casm.bat