md
set mem
800
md
vb off

; multireference alignment and classification
; uses the following procedures:
; alqr.mar
; centr.mar
; combat.mar



; number of images
x19=4467
; image size
x20=75
; file code
x22=6


; radius for alignment and mask
x21=35
; first radius for alignment
x26=5
; Search range 
; (the search for translation parameters will be 
; restricted to +/- search range).
x27=2
x28=1
; number of factors in CA S
x23=30
; number of groups
x24=100
; number of iterations of the multireference alignment
x25=5

; Prepare selection files for APMQ
DOC CREATE
mares/selectg
1
1-x24

DOC CREATE
mares/select
1
1-x19

;  Center all the images
;@centr[x19]
;../win/ser*****
;center/scr****
;mares/shiftr001

x91=0
do lb18 x92=1,x19
sd x92,x91,x91,x91,x91
mares/trns001
lb18
sd e
mares/trns001

; create mask for PCA
mo
_9
x20,x20
C
x21

; copy all the centered images to the stack
do lb21 x31=1,x19
cp
center/scr{****x31}
mares/data001@{****x31}
lb21
cp
mares/data001@
mares/data002@
;
; Iterate mutlireference alignment
do lb8 x81=1,x25
x82=x81+1
; Create SEQ file for PCA
ca si
mares/data002@****
1-x19
_9
SEQ{***x22}
Y
(0.0)
; Do the PCA
ca s
SEQ{***x22}
x23,x22
Y
0
; create subdirectory for the classification results
; and the averages
vm
echo "if(-d cluster{**x81}) rm -r cluster{**x81}"|csh
vm
mkdir cluster{**x81}
;
; classification
cl km
(2)
x22
x24
1-x23
(0)
(0)
cluster{**x81}/sel***
cluster{**x81}/cluster

; Get the averages
do lb1 x71=1,x24
as dc
mares/data002@****
cluster{**x81}/sel{***x71}
a
cluster{**x81}/avg{***x71}
cluster{**x81}/var{***x71}
lb1


; Optional
; Put averages in a decent orientation
@centr[x24]
cluster{**x81}/avg***
cluster{**x81}/cvg***
cluster{**x81}/centar


; Align averages roationaly
ap ra
cluster{**x81}/cvg***
1-x24
x26,x21
(0)
F
cluster{**x81}/haldoc

do lb41 x90=1,x24
ud ic,x90,x41,x42,x43,x44
cluster{**x81}/haldoc
rt sq
cluster{**x81}/cvg{***x41}
cluster{**x81}/rvg{***x41}
x42
(0,0)
lb41
ud ice
cluster{**x81}/haldoc

; end of optional alignment of the averages


; do the multireference alignment
ap mq
cluster{**x81}/rvg***
mares/selectg
x27
x26,x21
mares/data002@****
mares/select
cluster{**x81}/apmq001




@alqr
cluster{**x81}/apmq001
mares/data001@****
mares/data002@****
mares/trns{***x81}
mares/trns{***x82}

lb8

en