NB. A program to apply Saaty's method in an interactive fashion.
NB. according to the link in http://actifeld.com/A%20Possible%20Method.doc

NB. 2006-03-02
NB. John Randall developed the basic algorithm on the J Programming Forum
NB. Raul Miller, Tarmo Veskioja, Devon McCormick, and Roger Hui contributed to the dialog
NB. Bill Harris put a few other pieces on.

NB. The steps:
NB. Define the Client Preference Matrix as described on p. 6 of that article
NB. Run (%+/)power on that matrix to find the principle eigenvector.
NB. Run cr to calculate the Consistency Ratio; if it is much in excess of
NB.   0.1, the judgments used to create the CPM are probably too random.
NB. Create preference matrices for each of the alternatives on each of 
NB.   the dimensions used in the CPM.  Calculate the principal eigenvectors
NB.   and consistency ratio for each.
NB. Create the Option Performance Matrix by assembling all the principal 
NB.   eigenvectors for the preference matrices.
NB. Multiply OPM time CVV to get the Value For Money vector.  

NB. Calculate principal eigenvectors of the matrix cpm by
NB.   (%+/) power cpm

NB. Calculate the principal eigenvalue of the matrix cpm by
NB.   ev cpm

NB. Calculate the relative consistency ratio of the matrix cpm by
NB.   cr cpm 

NB. The VFM shows in each position the relative preference the decision-maker
NB.   has shown for that alternative: OPM * CVV = VFM

NB. If the Option Performance Matrix is
NB. opm=: |:a,b,c,:d  
NB. where a, b, c, and d are the principal eigenvectors of each option matrix

NB. then 
NB. (%/+)opm mp cpm 
NB. is the Value For Money (VFM) matrix, giving the relative 
NB. values of each alternative

NB. The power method for calculating the principal eigenvector
mp=:+/ . *
normalize=:%(>./@:|)
iterate=:normalize@:mp
init=:[: ? # # 0:
power=:13 : 'y&iterate^:_ init y'


NB. Calculate the principal eigenvalue
NB. Postmultiply the current principal eigenvector by the matrix
NB. Divide that vector (unnormalized) by the principal eigenvector
NB. See p. 258 of Kincaid and Cheney's Numerical Analysis 3rd edition
cureigenvec =: (%+/) @: power
neweigenvec =: mp cureigenvec
ev=: [: {. neweigenvec % cureigenvec

NB. Calculate the absolute consistency index
ci =: (ev - #) % <: @: #

NB. Calculate the consistency ratio
NB. MACI is the mean absolute consistency index of random matrices
am=: +/%#
maci =: 13 : 'am ci"2 st1"0 (x # y)'
NB. maci_precomputed =: 1000 maci"0 (3+i.8)   NB. Tarmo claims it's exactly 0.5245 for 3
index =: <:^:3 & #
maciv =: 13 : '(index y) { maci_precomputed' NB. Get the value of 
cr=: ci % maciv

NB. .........................................................................
NB. Random Saaty matrices for testing
NB. st is Tarmo Veskioja's verb (e.g., st 4)
NB. st1 is Roger Hui's faster but more space-intensive equivalent
saaty_scale =: \:~ }. (,%)~ >:i.9
str =: 13 : '( =/~ i.#y) + (y*bd ) + |: (bd=. >/~ i.#y) * %y'
saaty_rand =: 13 : '(y , y) $ saaty_scale {~ ? (#saaty_scale) #~ *: y'
st =: 13 : 'str saaty_rand y'        NB. st 4
st1=: 3 : 0
 i=. 1+(,~y) ?@$ #saaty_scale
 b=. >/~i.y
 ((b*i) + |:b*i+#saaty_scale) { 1,(,%) saaty_scale
)

NB. Test matrices from Geoff Coyle's paper referenced above:
NB. peigen s/b 0.232,0.402,0.061,0.305, cr of 0.055
cpm=: >(1,(%3),5, 1);(3,1,5,1);((%5),(%5),1,(%5));(1,1,5,1)

NB. peigen s/b 0.751,0.178,0.071 w/ cr of 0.072
exp=: >(1,5,9);((%5),1,3);((%9),(%3),1) NB. Expenses

NB. peigen s/b 0.480,0.406,0.114 w/ cr of 0.026
und=: >(1,1,5);(1,1,3);((%5),(%3),1) NB. "Understandability

NB. peigen s/b 0.077,0.231,0.692 w/ cr of 0
rod=: >(1,(%3),%9);(3,1,%3);9 3 1 NB. "Replication of detail", p. 7

NB. peigen s/b 0.066,0.615,0.319 w/ cr of 0
pod=: >(1,(%9),(%5));(9,1,2);(5,(%2),1) NB. Prediction of dynamics




scan =: 3 : 0
 for_j. i.#x do.
	sum =. sum + (x matrices"0 y) % (i.17^x)
 end.
)