#Example from: Gradient Projection Algorithms and Software for 
#  Arbitrary Rotation Criteria in Factor Analysis.
#  by Coen A. Bernaards and Robert I. Jennrich
#  Website: http://www.stat.ucla.edu/research

 Sys.getenv("R_LIBS")
 library()
 require("GPArotation")
 search()
 Sys.info()

 data("Thurstone", package="GPArotation") 
 if (!exists("box20")) stop("Test data not found. Testing stopped.")

fuzz <- 1e-5 
all.ok <- TRUE  

#  Thurstone's box problem.  (1947, p. 136) 
# The matrix box20 is the initial loading matrix from Thurstone's box problem.
            
# This takes a lot of iterations to converge at a higher tolerance
qbox20 <- quartimax(box20, eps=1e-5)
qbox20G <- GPForth(box20, Tmat=diag(1,3), method="quartimax", eps=1e-5)

 if( fuzz < max(abs(qbox20$loadings - qbox20G$loadings))) {
    cat("Calculated value is not the same as test value in test Thurstone 1. Value:\n")
    print(qbox20$loadings - qbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qbox20$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

#qbox20$Th - qbox20G$Th 

# These values compare with those in:
#    http://www.stat.ucla.edu/research/web.pdf
 tst <- t(matrix(c(
   0.0104916072210123716, -0.993396087928394733, -0.089861775335686706,
   0.1584646383898045685, -0.167305085570175344, -0.967087879524061056,
   0.9822741057703969769, -0.094961339079248266, -0.081938545344928893,
   0.1249962020162782989, -0.597065497283680413, -0.789290657131387352,
   0.8695614167874907707, -0.471622450093366785, -0.090438968384549553,
   0.8757114893176747294, -0.141012080768234127, -0.452333925937943637,
   0.0679423211019681700, -0.811411071716238719, -0.588554936857709099,
   0.4066768108416509708, -0.907862149146695163, -0.115673202040957226,
   0.5770808894249742638, -0.142370726163931066, -0.806527261406603468,
   0.1012712863762783577, -0.723336747696182614, -0.694640249329106285,
   0.5000928657774492692, -0.949746569049947253, -0.046846346456817907,
   0.7412589798326677526, -0.140350561965914555, -0.663578062154924320,
   0.0055655501003109590, -0.983847100401775698, -0.120037109608235590,
   0.2142330103903098415, -0.119429100752156334, -0.947421187831809397,
   0.9550804066106526324, -0.108275659756619305, -0.039227521113362487,
   0.7823218737697450464, -0.405437596810190704, -0.439275358874331168,
   0.3626971102221024923, -0.753122462957226402, -0.546281394544768872,
   0.0162483298780003657, -0.966230359337758582, -0.052114148464710915,
   0.1076692386876715729, -0.206734953950642314, -0.934620775424686911,
   0.9744239420161749932, -0.092650552854598708, -0.090828719474599584
   ), 3, 20))
 
 if( fuzz < max(abs(qbox20$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 2. Value:\n")
    print(qbox20$loadings, digits=18)
    cat("difference:\n")
    print(qbox20$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

 tst <- t(matrix(c(
   0.57232345894276127, -0.60751194947821441, -0.55079496147384377,
   0.60249460283341838,  0.76716797198365361, -0.22012168525406509,
   0.55627880770383020, -0.20587018726291534,  0.80509089803322043
   ), 3, 3))
 
 if( fuzz < max(abs(qbox20$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 3. Value:\n")
    print(qbox20$Th, digits=18)
    cat("difference:\n")
    print(qbox20$Th - tst, digits=18)
    all.ok <- FALSE  
    } 


# sorted absolute loading plots. 
sal <- abs(c(loadings(qbox20)))[order(abs(c(loadings(qbox20))))]
plot(seq(length(sal)), sal)


 #compare quartimax rotation of the initial loading matrix box20.

 if( fuzz < max(abs(loadings(qbox20) - box20 %*% qbox20$Th ))) {
    cat("Calculated value is not the same as test value in test Thurstone 4. Value:\n")
    print(loadings(qbox20), digits=18)
    cat("difference:\n")
    print(loadings(qbox20) - box20 %*% qbox20$Th, digits=18)
    all.ok <- FALSE  
    } 


            
qminbox20G <- GPFoblq(box20, Tmat=diag(1,3), method="quartimin", eps=1e-5)
qminbox20  <- quartimin(box20, eps=1e-5)

 if( fuzz < max(abs(loadings(qminbox20) - qminbox20G$loadings))) {
    cat("Calculated value is not the same as test value in test Thurstone 5. Value:\n")
    print(qminbox20G$loadings , digits=18)
    cat("difference:\n")
    print(loadings(qminbox20) - qminbox20G$loadings, digits=18)
    all.ok <- FALSE  
    } 

#qminbox20$Th - quartimin(box20)$Th 

# These values compare with those in:
#    http://www.stat.ucla.edu/research/web.pdf
 tst <- t(matrix(c(
   -0.099561899210599963, -1.0236437309424475384,  0.017110338313848200,
   -0.007103778102200991,  0.0427848301281630802, -1.009962780073245581,
    1.012864497258948226,  0.0331727792925069487,  0.050367710973030555,
   -0.054843850612513692, -0.4493155290974688021, -0.772334543778026350,
    0.856287122381722998, -0.3740197232441037078,  0.069350368268248391,
    0.835580575619599641,  0.0487450425576793633, -0.360381644212301344,
   -0.102893671454670210, -0.7226715938020771279, -0.537456650126404090,
    0.322103633211960838, -0.8816846447967544576,  0.031159743715387874,
    0.462799683447739529,  0.0852338438217692257, -0.783762970578423479,
   -0.076585435689138226, -0.6043060025891554554, -0.658295846696152820,
    0.427772530893690217, -0.9288687512327726825,  0.122866182561916254,
    0.659408232467282085,  0.0772080094990600374, -0.607348040513722709,
   -0.108761719100651882, -1.0079608432113262850, -0.017378089000366713,
    0.059518597564186392,  0.0955950614351480238, -0.986779686330629513,
    0.989890866913205381,  0.0071520817823045348,  0.094691644950703049,
    0.713733277219835149, -0.2427293600063723522, -0.328268187306521242,
    0.220344503737931546, -0.6353746612195683152, -0.459661643730432223,
   -0.084703580704062989, -1.0022284232457450148,  0.055740317456252478,
   -0.059151779416785115, -0.0113377397453605679, -0.976867596293413132,
    1.003360458549731771,  0.0365098037316876067,  0.039427150580815938
    ), 3, 20))
 
 if( fuzz < max(abs(qminbox20G$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 6. Value:\n")
    print(qminbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qminbox20G$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 

 tst <- t(matrix(c(
   1.00000000000000000, -0.25676300454795098, -0.32155119431295237,
  -0.25676300454795098,  1.00000000000000000,  0.33656790396842257,
  -0.32155119431295237,  0.33656790396842257,  1.00000000000000000
   ), 3, 3))
 
 if( fuzz < max(abs(qminbox20G$Phi - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 7. Value:\n")
    print(qminbox20G$Phi, digits=18)
    cat("difference:\n")
    print(qminbox20G$Phi - tst, digits=18)
    all.ok <- FALSE  
    } 


#To fuzz precision the rotated loading matrix and the factor cor-
#relation matrix phi are identical to those produced using the oblique GP
#algorithm with exact derivatives.

 if( fuzz < max(abs(qminbox20G$Phi -
                  t(qminbox20G$Th )%*% qminbox20G$Th ))) {
    cat("Calculated value is not the same as test value in test Thurstone 8. Value:\n")
    print(qminbox20G$Phi, digits=18)
    cat("difference:\n")
    print(qminbox20G$Phi - 
                  t(qminbox20G$Th )%*% qminbox20G$Th, digits=18)
    all.ok <- FALSE  
    } 

 #compare quartimin rotation of the initial loading matrix box20.
   if( fuzz < max(abs(qminbox20G$loadings - box20 %*% solve(t(qminbox20G$Th))))) {
    cat("Calculated value is not the same as test value in test Thurstone 9. Value:\n")
    print(qminbox20G$loadings, digits=18)
    cat("difference:\n")
    print(qminbox20G$loadings - box20 %*% solve(t(qminbox20G$Th)), digits=18)
    all.ok <- FALSE  
    } 

 data("box26", package="GPArotation") 
 if (!exists("box26")) stop("Test data box26 not found. Testing stopped.")

qbox26  <- GPForth(box26, Tmat=diag(1,3), method="quartimax", eps=1e-5)

 tst <- t(matrix(c(
   0.6245197355925140581, -0.2708954695931116152,  0.7151983951389878635,
   0.7386116884036847408,  0.6266342260884526505, -0.0617439911892987553,
   0.7803093788467402314, -0.3830982859243221017, -0.4578886072022986253,
   0.8540550453155928423,  0.2886436985992582027,  0.4062915145925659610,
   0.8810593765418006651, -0.4428658074662961130,  0.1233946983666596581,
   0.9084731768740617053,  0.1540526132602804965, -0.3723026715563940159,
   0.8150592858039771293,  0.0441965358534676597,  0.5600768044145943980,
   0.8466584455973064083,  0.4551177395514792168,  0.1889929089788950356,
   0.8156808837280125069, -0.4090629943132625956,  0.3690652552112651530,
   0.9629492340906220527, -0.4781483041690369196, -0.0866081507974762743,
   0.8731366884896356595,  0.3451069860590937899, -0.2914969834947889749,
   0.8921854600753849063, -0.0276323108621970258, -0.4257376659710629951,
  -0.0938760381595044741, -0.7873218033841372643,  0.6012450975895150540,
   0.0938760381595044741,  0.7873218033841372643, -0.6012450975895150540,
  -0.0986092863860908303,  0.1513605567468480073,  0.9692559984337008050,
   0.0986092863860908303, -0.1513605567468480073, -0.9692559984337008050,
  -0.0189573629854957251,  0.9527983290277913797,  0.2944078167958268377,
   0.0189573629854957251, -0.9527983290277913797, -0.2944078167958268377,
   0.8394181189595459891,  0.3631767908642606346,  0.3398717995655929913,
   0.8703065201362156778, -0.4691145408161159214,  0.0770980453920554615,
   0.9141063746617547059,  0.1583184861345137973, -0.3535658252020681958,
   0.8348118627305495254,  0.3535663452183119837,  0.3271666140872500073,
   0.8541352373790773722, -0.4476738735312740247,  0.0569042988261160704,
   0.9034738474019414767,  0.1663655738425987851, -0.3227406124130587362,
   0.9861758757457432800,  0.0103496363116840455,  0.0635926656567585569,
   0.9643516568468981642,  0.0660181478622221818, -0.0304218028637989850
   ), 3, 26))
 
 if( fuzz < max(abs(qbox26$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 10. Value:\n")
    print(qbox26$loadings, digits=18)
    cat("difference:\n")
    print(qbox26$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
   0.9996572020207266096, 0.0216275672176080257,  0.0147555679097727491,
  -0.0158190757965277796, 0.9480178905874908635, -0.3178235925273457108,
  -0.0208622934749700742, 0.3174812237948764770,  0.9480350400953921897
   ), 3, 3))
 
 if( fuzz < max(abs(qbox26$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 11. Value:\n")
    print(qbox26$Th, digits=18)
    cat("difference:\n")
    print(qbox26$Th - tst, digits=18)
    all.ok <- FALSE  
    } 


qminbox26 <- GPFoblq(box26, Tmat=diag(1,3), method="quartimin", eps=1e-5)

 tst <- t(matrix(c(
   0.6088436426802223966, -0.2567107018725688361,  0.7213648290819488773,
   0.7318447535507376367,  0.6298398026581654152, -0.0549983771960348838,
   0.7973321695017724364, -0.3855960314746548212, -0.4504478973568259437,
   0.8392144987741166906,  0.2994932968625432235,  0.4143558581243267924,
   0.8833452352200144020, -0.4361046712803113290,  0.1319331147095905710,
   0.9161366872228343672,  0.1535557844336666034, -0.3638337328539109072,
   0.7993355454002614158,  0.0571270784641514026,  0.5678963531379384033,
   0.8354288250614068101,  0.4626764152757318893,  0.1968789749765105790,
   0.8109923806202916641, -0.3989909333845649830,  0.3770226870580207779,
   0.9712737747877250305, -0.4740722765307348041, -0.0773243882106463137,
   0.8761501947960563808,  0.3456235893514668089, -0.2834183138879167174,
   0.9036601763684347643, -0.0290211959776035672, -0.4173652812159966974,
  -0.0995525797764766768, -0.7788574612781464790,  0.6007791331268093060,
   0.0995525797764766768,  0.7788574612781464790, -0.6007791331268093060,
  -0.1264036712449473909,  0.1653130238928011975,  0.9684661160120416890,
   0.1264036712449473909, -0.1653130238928011975, -0.9684661160120416890,
  -0.0392946742598458687,  0.9571059478962877787,  0.2939285303852590125,
   0.0392946742598458687, -0.9571059478962877787, -0.2939285303852590125,
   0.8253744379910458173,  0.3729516010405902748,  0.3477554718030251846,
   0.8741734789142978634, -0.4631063486451737488,  0.0855349365396926159,
   0.9212130243051569467,  0.1581334796580046720, -0.3450412531516501846,
   0.8212340853954427367,  0.3631252613622908965,  0.3350076577679809153,
   0.8582635618771776720, -0.4420579024138228674,  0.0651757040961165046,
   0.9096561314838297330,  0.1665824736284239604, -0.3143133889989875307,
   0.9840845767693481294,  0.0168070160966761091,  0.0729425956763933708,
   0.9640420478016114014,  0.0709475796833391181, -0.0213192081807395371
   ), 3, 26))
 
 if( fuzz < max(abs(qminbox26$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 12. Value:\n")
    print(qminbox26$loadings, digits=18)
    cat("difference:\n")
    print(qminbox26$loadings - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
   1.000000000000000	 ,  0.00767934084449363279,  0.0170654511973979163,
   0.00767934084449363279,  1.000000000000000	  , -0.0144994900961642244,
   0.01706545119739791630, -0.01449949009616422445,  1.000000000000000     
    ), 3, 3))
 
 if( fuzz < max(abs(qminbox26$Phi - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 13. Value:\n")
    print(qminbox26$Phi, digits=18)
    cat("difference:\n")
    print(qminbox26$Phi - tst, digits=18)
    all.ok <- FALSE  
    } 



 tst <- t(matrix(c(
    0.9993401424148040668, 0.0347479564402226465,  0.0408645923859655008,
   -0.0179660947915933414, 0.9476477730670300748, -0.3324117322929439067,
   -0.0315673755054017430, 0.3174212937474846785,  0.9422486536594960604
    ), 3, 3))
 
 if( fuzz < max(abs(qminbox26$Th - tst ))) {
    cat("Calculated value is not the same as test value in test Thurstone 14. Value:\n")
    print(qminbox26$Th, digits=18)
    cat("difference:\n")
    print(qminbox26$Th - tst, digits=18)
    all.ok <- FALSE  
    } 



cat("tests completed.\n")

if (! all.ok) stop("some tests FAILED")


cat("tests completed.\n")

if (! all.ok) stop("some tests FAILED")