From de4cdf7e22ced8df46496f131ed2eab57b1246bd Mon Sep 17 00:00:00 2001 From: Frank E Harrell Jr Date: Mon, 28 Oct 2024 05:10:05 +0000 Subject: [PATCH] version 5.2-0 --- DESCRIPTION | 12 ++-- MD5 | 42 ++++++------ NAMESPACE | 2 +- NEWS | 9 ++- R/describe.s | 20 +++--- R/pMedian.r | 27 ++++++++ R/qcrypt.r | 3 + inst/tests/pMedian.r | 94 +++++++++++++++++++++++++++ man/Hmisc-internal.Rd | 1 + man/describe.Rd | 6 +- man/latex.Rd | 4 +- man/pMedian.Rd | 33 ++++++++++ man/popower.Rd | 6 +- man/transcan.Rd | 2 +- src/Hmisc.h | 8 +++ src/cidxcn.f | 88 ------------------------- src/cidxcn.f90 | 95 +++++++++++++++++++++++++++ src/cidxcp.f | 92 -------------------------- src/cidxcp.f90 | 94 +++++++++++++++++++++++++++ src/hlqest.f90 | 144 +++++++++++++++++++++++++++++++++++++++++ src/hoeffd.f | 120 ---------------------------------- src/hoeffd.f90 | 146 ++++++++++++++++++++++++++++++++++++++++++ src/init.c | 2 + src/jacklins.f | 23 ------- src/jacklins.f90 | 31 +++++++++ src/maxempr.f | 68 -------------------- src/maxempr.f90 | 81 +++++++++++++++++++++++ src/rcorr.f | 106 ------------------------------ src/rcorr.f90 | 123 +++++++++++++++++++++++++++++++++++ src/wclosest.f | 57 ----------------- src/wclosest.f90 | 72 +++++++++++++++++++++ 31 files changed, 1013 insertions(+), 598 deletions(-) create mode 100644 R/pMedian.r create mode 100644 inst/tests/pMedian.r create mode 100644 man/pMedian.Rd delete mode 100644 src/cidxcn.f create mode 100644 src/cidxcn.f90 delete mode 100644 src/cidxcp.f create mode 100644 src/cidxcp.f90 create mode 100644 src/hlqest.f90 delete mode 100644 src/hoeffd.f create mode 100644 src/hoeffd.f90 delete mode 100644 src/jacklins.f create mode 100644 src/jacklins.f90 delete mode 100644 src/maxempr.f create mode 100644 src/maxempr.f90 delete mode 100644 src/rcorr.f create mode 100644 src/rcorr.f90 delete mode 100644 src/wclosest.f create mode 100644 src/wclosest.f90 diff --git a/DESCRIPTION b/DESCRIPTION index da71c01..17fbc62 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: Hmisc -Version: 5.1-3 -Date: 2024-05-18 +Version: 5.2-0 +Date: 2024-10-26 Title: Harrell Miscellaneous Authors@R: c(person(given = "Frank E", @@ -14,7 +14,7 @@ Authors@R: email = "charles.dupont@vumc.org", comment = "contributed several functions and maintains latex functions")) Maintainer: Frank E Harrell Jr -Depends: R (>= 4.1.0) +Depends: R (>= 4.2.0) Imports: methods, ggplot2, cluster, rpart, nnet, foreign, gtable, grid, gridExtra, data.table, htmlTable (>= 1.11.0), viridis, htmltools, base64enc, colorspace, rmarkdown, knitr, Formula @@ -33,11 +33,11 @@ License: GPL (>= 2) LazyLoad: Yes URL: https://hbiostat.org/R/Hmisc/ Encoding: UTF-8 -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 NeedsCompilation: yes -Packaged: 2024-05-18 13:26:32 UTC; harrelfe +Packaged: 2024-10-27 14:33:00 UTC; harrelfe Author: Frank E Harrell Jr [aut, cre] (), Charles Dupont [ctb] (contributed several functions and maintains latex functions) Repository: CRAN -Date/Publication: 2024-05-28 07:10:02 UTC +Date/Publication: 2024-10-28 05:10:05 UTC diff --git a/MD5 b/MD5 index 934d1ed..58b6816 100644 --- a/MD5 +++ b/MD5 @@ -1,7 +1,7 @@ 16f8ea43efba90112351c5ecb33021c4 *COPYING -77dbeb98ec3ce4c4f0f2afdd5144db1a *DESCRIPTION -567cedef1920aa3944bac88ae5e528ff *NAMESPACE -39f3a27fdffa3badf4b28934948e4064 *NEWS +16e2d134d31f23bb056c142bd79ed34d *DESCRIPTION +5155c825bd5b38425b18e7928c7dfefb *NAMESPACE +eb8c726e113ab88278b6f2198926acc5 *NEWS 76f90acbe9bf8ad2eaaa6b8813eb7c82 *R/AFirst.lib.s 8fdbfa3b900cbbb42d94898093ec7fa5 *R/Cs.s 314b24f935cf8e2ee7a677ed26479b83 *R/GiniMd.s @@ -35,7 +35,7 @@ c04e31869fdc5107eb6af6886eadb566 *R/cpower.s 97dc5dcc48cb719088b817833067352c *R/dataRep.s ec8af558be91e1fa3415f1a99dd469b2 *R/dates.s 3f02d2486d14e095a96fe5ee10e483c7 *R/deff.s -1356fda7b10bcc8bb85c32b64e25896f *R/describe.s +0cc7d70f97385be3c9e5942357e4b4e0 *R/describe.s 6cb5b3a77860b605f102bb528e84a071 *R/discrete.s 3c7e7885dccb0e283f96f91fe1b170ba *R/dotchart3.s 3d1c3094a57c0fcf886f87367b539d86 *R/dotchartpl.s @@ -102,6 +102,7 @@ dd4806066d614756cd8be2bef6bad3dd *R/na.keep.s abc0f1c38d72a9bf3e59ce9d194c71ee *R/nobsY.s 27019e9e90730ac0b25014b51d68ec66 *R/nstr.s cf2866606281d01f39e99d3205f50ae3 *R/num.intercepts.s +6dedd1730fffe3aa6d68658a227fabf0 *R/pMedian.r a50c8048a7546e224ca402d362c084cd *R/pairUpDiff.r 215f9a06ffc6df02fb67576ad0da83b9 *R/panel.abwplot.s 68d8bbb7b9e6f5fa975b37f2e1245291 *R/panel.bpplot.s @@ -115,7 +116,7 @@ a9771408bed89b32b030653ea46a614f *R/princmp.r 9a1119c7949407c39968645e921a87a6 *R/print.char.list.s 5e6df449d8b6cd6f06693d7298d47150 *R/printL.r efe9a9a6c90594d1e7425cc46f2fc061 *R/pstamp.s -21e5f0670676be1ba7b4b2ad135de5e7 *R/qcrypt.r +e2cf58cc149b7d657e44a976f5e8f185 *R/qcrypt.r f3c43a8d39f1458c92a8d2d7df290d9a *R/r2describe.r 2679a3e28310299f523e0ffb440db739 *R/rcorr.cens.s b9759f7459856f6eeb85e4ee3832f7ad *R/rcorr.s @@ -225,6 +226,7 @@ a40f892874e68edf0c621971b200260d *inst/tests/latexTabular.r b0fc6d4900270c13a20eafd31293bc9c *inst/tests/latexpng.r 0c6b7314112171894df7511e12398f05 *inst/tests/mChoice.r c0668a4323dcc6d83736668c00456d23 *inst/tests/minor.tick.r +de29fe314e56e44256137d61652066e4 *inst/tests/pMedian.r ca3aeb15387cea45fb60b35bdd1b572b *inst/tests/panelbp.r 39662c6cbeab67ca1544cb736ea64f70 *inst/tests/plot.summaryM.plotly.r badcf6388bfae0071da3c10a2b77bac6 *inst/tests/procmeans.txt @@ -255,7 +257,7 @@ e1646179a6487ec47f043ba1125a5438 *inst/tests/wtd.r 7e9e0d3131e1e1fc0e423359e1c829f6 *man/Cs.Rd cec2ab865d8b920cb20e4157dc2a3faf *man/Ecdf.Rd ac4fff37fde2e6d0feae41ceda669399 *man/GiniMd.Rd -2b084c69d1db718784be029130560d30 *man/Hmisc-internal.Rd +a40b9629c493e0d79c8d615968b05b45 *man/Hmisc-internal.Rd a8d3213ffdeb74c3f50bb2d54a8b9ab5 *man/Lag.Rd 352fb4122854747013fc8d3c96ddb360 *man/Merge.Rd 5924598d446057203533bc6ee2970150 *man/Misc.Rd @@ -290,7 +292,7 @@ f76ccfc1d3e0558d9be0a4f89d6f8f62 *man/cut2.Rd 3eb4090f0c25a33dcf3452294ca80ec1 *man/data.frame.create.modify.check.Rd 30aa19bbc67a730bfe65d71ecd473bb4 *man/dataRep.Rd 1df7d66a5fcdc97254d4e70a8d487af4 *man/deff.Rd -6edbbe5d2ba26d80161640c46971f853 *man/describe.Rd +6ec474869dfc2350b2d6b178751f8afc *man/describe.Rd 4c2942aab548faa49c5a389117b48b70 *man/discrete.Rd 94748a0c9bb7d0290a8b2dc7f4a0bb7d *man/dotchart2.Rd 67f28b3925d62e99e612e36532d3149f *man/dotchart3.Rd @@ -338,7 +340,7 @@ b785dbffdb86727e8650cf7e19df0377 *man/intMarkovOrd.Rd 4831b433c139db936cecf4dff76b9a05 *man/labcurve.Rd fd6547343482e782269052b8ed6cc1c5 *man/label.Rd d0d1b03b264e6d6456c4e03e5e793c3f *man/latestFile.Rd -bebcb911d2f805ec0079cf0d7e4c02c7 *man/latex.Rd +99bcdc4206da2fedc3111efc1310bd97 *man/latex.Rd 9c1dfa5242426105b1ca4669b0977997 *man/latexCheckOptions.Rd c4a248262270a156d9723212daba90ff *man/latexDotchart.Rd 32f49f9726f15c0cae64179b3df74f20 *man/latexTabular.Rd @@ -363,6 +365,7 @@ f9eb6868ac6e099103d654c1c1084b21 *man/nin.Rd be66aaec2a64401335915c4671d50918 *man/nobsY.Rd 2f7e82883fd0be5a86a3bec4d4dc250a *man/nstr.Rd 1d560c844c38c8f7289bb27e21efe171 *man/num.intercepts.Rd +5aed2c0b67a2c5c52af1d39e60bcdb94 *man/pMedian.Rd 76b1975307ab89cc2eb19f12040342ce *man/pairUpDiff.Rd ab8edbb03663977e40ac2ffa3025d39f *man/panel.bpplot.Rd b8247448b0c2c6d876ce693f23c9fe85 *man/partition.Rd @@ -372,7 +375,7 @@ fd3bfeefd15a79c301f02a7c4631f8d1 *man/plotCorrM.Rd a3aed45d350c39121d61700f94e859cb *man/plotCorrPrecision.Rd 88d93c5117d432cbe4ba996c5c74d0d8 *man/plotlyM.Rd 81937c2b724d0ec42f817ebb88748e70 *man/plsmo.Rd -66c7af813ac11f2bd92aebcc6ff33d3e *man/popower.Rd +6a480cdaeb1dbf34525224011cbe17b9 *man/popower.Rd ddba6cfbe8781522fad083b7bc2c8b92 *man/princmp.Rd 52be83076b1f0442cb5ab5c15bb16a3f *man/print.char.list.Rd 0a25e1c6f349789159b7239ae0f0f620 *man/print.char.matrix.Rd @@ -435,7 +438,7 @@ d9e66a07144bec4e7365da937119b5ea *man/symbol.freq.Rd 996241a9c13b901a3d5133de1624b9ed *man/testCharDateTime.Rd ec1a45cd81fe0bb55b45eb5b0ead1e73 *man/tex.Rd 6a55f6a9b70721cec9c3b8830bceaa0e *man/transace.Rd -b598e89a3f72832bbe5a0e9dafb7c7c8 *man/transcan.Rd +befd3ceb2cd8cec372ab9c079cbbac2a *man/transcan.Rd 7e9eae90e12d52f637f46ce78bfc0c53 *man/translate.Rd 76093ee37d83f9e13c9be04fbaaca9bd *man/trunc.POSIXt.Rd 2a6ef4d194b66158058fcfe0055d718f *man/units.Rd @@ -453,15 +456,16 @@ a7ee5852bd965cedb06b50e051df63b1 *man/windows/sas.get.Rd c63b3fa2193829f494e39d8cd557cd83 *man/yearDays.Rd c87ff598535f37eaf4c78f044e722b74 *man/ynbind.Rd 8360e406703130bb09302a64a9638860 *src/Hmisc.c -a8395ba413d746154d7e89ce07e1e11b *src/Hmisc.h -06d9b6e7e6a7e22898fc0134691889b3 *src/cidxcn.f -ac182f4144cfd7853f252be1e4b3d75b *src/cidxcp.f -08ccfbf5a4efcc18c629a8f91abc6b78 *src/hoeffd.f -781a196766abf3ebddc10bc658ff22a2 *src/init.c -d6bc8b5ebd92a27f1c9e3522d46fa466 *src/jacklins.f +feacb897c62776acb7e8ebffa7446262 *src/Hmisc.h +72742560872c7a863b1aa9ce21160d4a *src/cidxcn.f90 +a519b7e05097490440cb744c57392de4 *src/cidxcp.f90 +1d61e62a27b2830ad58ac6645b11c5df *src/hlqest.f90 +3e94078aa23d484e0bd94375b813855e *src/hoeffd.f90 +23773dc91efdd1871530587df857ff02 *src/init.c +1f84fad7c8d73e5839030ce4d2cf1f85 *src/jacklins.f90 10538b0cf98c499976a5bb6d74a5195e *src/largrec.f 2d0aa4a90313d361f628383dd5fcf7fb *src/mChoice.c -703f73252fe6efaf8244b0713aaabe6b *src/maxempr.f +9e3c72bd4e48c78170545e5e0c50d7be *src/maxempr.f90 69d3a08b2cd911161ba22af68f97f611 *src/nstr.c 5fd3cecbf580b3cb32f365f1f9036794 *src/ranksort.c df53032c14b9c8a41bfe9529411c7729 *src/ratfor/cidxcn.r @@ -471,7 +475,7 @@ eec5235344d56eaaaff52805e71fdf30 *src/ratfor/jacklins.r daae316a69a5b013bdf334c4302725cd *src/ratfor/maxempr.r 29402d28cea92baf7631bc7ba15be2d2 *src/ratfor/rcorr.r 1e330cf030042a8aef488a254a6b1074 *src/ratfor/wclosest.r -def0335461e456fdd3d966856a7a781c *src/rcorr.f +1be4ba1903bae2d80ae7fb0b7e97a0b3 *src/rcorr.f90 32839fee8ce62db43ef9b527256a8737 *src/sas/exportlib.sas 6391fc48464e9b488b3626eaacbe967d *src/string_box.c -46912e086e55a7d273443ff5a6e37bde *src/wclosest.f +c6c9ebc7448e3c3fafdbac4c610344d4 *src/wclosest.f90 diff --git a/NAMESPACE b/NAMESPACE index ae44330..25accae 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -export("%nin%",.q,abs.error.pred,addggLayers,addMarginal,all.digits,all.is.numeric,approxExtrap,areg,areg.boot,aregImpute,aregTran,arrGrob,as.discrete,asNumericMatrix,ballocation,bezier,binconf,biVar,bootkm,bpower,bpower.sim,bpplot,bpplotM,bpplt,bppltp,bpx,bsamsize,bystats,bystats2,capitalize,catTestchisq,Cbind,ceil,character.table,chiSquare,ciapower,cleanup.import,clowess,cnvrt.coords,code.levels,colorFacet,combine.levels,combineLabels,completer,combplotp,confbar,consolidate,contents,conTestkw,cpower,Cs,csv.get,cumcategory,curveRep,curveSmooth,cut2,datadensity,dataDensityString,dataframeReduce,dataRep,ddmmmyy,deff,describe,discrete,dhistboxp,dotchart2,dotchart3,dotchartp,dotchartpl,Dotplot,drawPlot,dvi,dvigv,dvips,ebpcomp,ecdfpM,Ecdf,ecdfSteps,equalBins,errbar,escapeBS,escapeRegex,estSeqMarkovOrd,estSeqSim,event.chart,event.convert,event.history,expr.tree,extractlabs,fillin,fImport,find.matches,first.word,fit.mult.impute,format.df,format.pval,formatCats,formatCons,formatDateTime,formatdescribeSingle,formatSep,formatTestStats,ftupwr,ftuss,Function,gbayes,gbayes1PowerNP,gbayes2,gbayesMixPost,gbayesMixPowerNP,gbayesMixPredNoData,gbayesSeqSim,get2rowHeads,getabd,getHdata,getLatestSource,GetModelFrame,getRs,getZip,ggfreqScatter,ggplotlyr,GiniMd,Gompertz2,groupn,grType,hashCheck,hdquantile,hidingTOC,hist.data.frame,histbackback,histboxp,histboxpM,histSpike,histSpikeg,hlab,hlabs,hoeffd,html,htmlGreek,htmlSN,htmlSpecial,htmlSpecialType,htmltabv,htmlTranslate,htmlVerbatim,importConvertDateTime,improveProb,impute,impute.transcan,inmChoice,inmChoicelike,intMarkovOrd,inverseFunction,invertTabulated,is.discrete,is.imputed,is.mChoice,is.present,is.special.miss,james.stein,jitter2,keepHattrib,Key,Key2,km.quick,knitrSet,labcurve,label,Label,labelLatex,labelPlotmath,Lag,largest.empty,latestFile,latex,latexBuild,latexCheckOptions,latexDotchart,latexNeedle,latexSN,latexTabular,latexTherm,latexTranslate,latexVerbatim,list.tree,llist,lm.fit.qr.bare,Load,Lognorm2,logrank,lookupSASContents,lrcum,makeNames,makeNstr,makeSteps,mApply,markupSpecs,mask,match.mChoice,matchCases,matrix2dataFrame,matxv,mbarclPanel,mbarclpl,mChoice,mdb.get,Mean,medvPanel,medvpl,meltData,Merge,mgp.axis,mgp.axis.labels,mhgr,minor.tick,monotone,monthDays,movStats,mtitle,multEventChart,multLines,na.delete,na.detail.response,na.include,na.keep,na.pattern,na.retain,naclus,nafitted.delete,Names2names,naplot,napredict.delete,naprint.delete,naprint.keep,naresid.delete,naresid.keep,nCoincident,nFm,nmChoice,nobsY,nomiss,nstr,num.denom.setup,num.intercepts,numeric.string,numericScale,oPar,optionsCmds,ordGridFun,ordTestpo,outerText,pairUpDiff,panel.bpplot,panel.Dotplot,panel.Ecdf,panel.plsmo,panel.xYplot,parGrid,partition.matrix,partition.vector,pasteFit,pBlock,pc1,plotCorrM,plotCorrPrecision,plotp,plotlyM,plotlyParm,plotlySave,plotmathTranslate,plotMultSim,plotpsummaryM,plsmo,pngNeedle,pomodm,popower,posamsize,prepanel.Dotplot,prepanel.Ecdf,prepanel.xYplot,princmp,printL,print.char.matrix,prList,prn,propsPO,propsTrans,prselect,prType,pstamp,putHcap,putHfig,putKey,putKeyEmpty,qcrypt,Quantile,Quantile2,r2describe,R2Measures,rcorr,rcorr.cens,rcorrcens,rcorrp.cens,rcspline.eval,rcspline.plot,rcspline.restate,rcsplineFunction,read.xportDataload,readSAScsv,redun,reformM,reLabelled,rendHTML,replace.substring.wild,reShape,responseSummary,restoreHattrib,rlegend,rlegendg,rm.boot,rMultinom,roundN,roundPOSIXt,runifChanged,runParallel,samplesize.bin,sas.codes,sas.get,sas.get.macro,sasdsLabels,sasxport.get,Save,scat1d,score.binary,sedit,sepUnitsTrans,seqFreq,setParNro,setTrellis,show.col,show.pch,showPsfrag,simMarkovOrd,simplifyDims,simPOcuts,simRegOrd,sKey,smean.cl.boot,smean.cl.normal,smean.sd,smean.sdl,smearingEst,smedian.hilow,solvet,somers2,soprobMarkovOrd,soprobMarkovOrdm,spearman,spearman.test,spearman2,spikecomp,spower,spss.get,src,stat_plsmo,stata.get,StatPlsmo,stepfun.eval,stratify,strgraphwrap,string.bounding.box,string.break.line,stringDims,stripChart,subplot,substi,substi.source,substring.location,substring2,summarize,summaryD,summaryDp,summaryM,summaryP,summaryRc,summaryS,symbol.freq,sys,t.test.cluster,table_formatpct,table_freq,table_latexdefs,table_N,table_pc,table_trio,tabulr,termsDrop,testCharDateTime,testDateTime,tex,tobase64image,transace,transcan,translate,trap.rule,trellis.strip.blank,truncPOSIXt,uncbind,units,unPaste,upData,upFirst,valueLabel,valueName,valueTags,valueUnit,var.inner,varclus,vlab,Weibull2,whichClosek,whichClosePW,whichClosest,wtd.Ecdf,wtd.loess.noiter,wtd.mean,wtd.quantile,wtd.rank,wtd.table,wtd.var,xInch,xless,xy.group,xYplot,xySortNoDupNoNA,yearDays,yInch,ynbind,zoom,"[<-.discrete","consolidate<-","is.na<-.discrete","label<-","length<-.discrete","substring2<-","valueLabel<-","valueName<-","valueTags<-","valueUnit<-") +export("%nin%",.q,abs.error.pred,addggLayers,addMarginal,all.digits,all.is.numeric,approxExtrap,areg,areg.boot,aregImpute,aregTran,arrGrob,as.discrete,asNumericMatrix,ballocation,bezier,binconf,biVar,bootkm,bpower,bpower.sim,bpplot,bpplotM,bpplt,bppltp,bpx,bsamsize,bystats,bystats2,capitalize,catTestchisq,Cbind,ceil,character.table,chiSquare,ciapower,cleanup.import,clowess,cnvrt.coords,code.levels,colorFacet,combine.levels,combineLabels,completer,combplotp,confbar,consolidate,contents,conTestkw,cpower,Cs,csv.get,cumcategory,curveRep,curveSmooth,cut2,datadensity,dataDensityString,dataframeReduce,dataRep,ddmmmyy,deff,describe,discrete,dhistboxp,dotchart2,dotchart3,dotchartp,dotchartpl,Dotplot,drawPlot,dvi,dvigv,dvips,ebpcomp,ecdfpM,Ecdf,ecdfSteps,equalBins,errbar,escapeBS,escapeRegex,estSeqMarkovOrd,estSeqSim,event.chart,event.convert,event.history,expr.tree,extractlabs,fillin,fImport,find.matches,first.word,fit.mult.impute,format.df,format.pval,formatCats,formatCons,formatDateTime,formatdescribeSingle,formatSep,formatTestStats,ftupwr,ftuss,Function,gbayes,gbayes1PowerNP,gbayes2,gbayesMixPost,gbayesMixPowerNP,gbayesMixPredNoData,gbayesSeqSim,get2rowHeads,getabd,getHdata,getLatestSource,GetModelFrame,getRs,getZip,ggfreqScatter,ggplotlyr,GiniMd,Gompertz2,groupn,grType,hashCheck,hdquantile,hidingTOC,hist.data.frame,histbackback,histboxp,histboxpM,histSpike,histSpikeg,hlab,hlabs,hoeffd,html,htmlGreek,htmlSN,htmlSpecial,htmlSpecialType,htmltabv,htmlTranslate,htmlVerbatim,importConvertDateTime,improveProb,impute,impute.transcan,inmChoice,inmChoicelike,intMarkovOrd,inverseFunction,invertTabulated,is.discrete,is.imputed,is.mChoice,is.present,is.special.miss,james.stein,jitter2,keepHattrib,Key,Key2,km.quick,knitrSet,labcurve,label,Label,labelLatex,labelPlotmath,Lag,largest.empty,latestFile,latex,latexBuild,latexCheckOptions,latexDotchart,latexNeedle,latexSN,latexTabular,latexTherm,latexTranslate,latexVerbatim,list.tree,llist,lm.fit.qr.bare,Load,Lognorm2,logrank,lookupSASContents,lrcum,makeNames,makeNstr,makeSteps,mApply,markupSpecs,mask,match.mChoice,matchCases,matrix2dataFrame,matxv,mbarclPanel,mbarclpl,mChoice,mdb.get,Mean,medvPanel,medvpl,meltData,Merge,mgp.axis,mgp.axis.labels,mhgr,minor.tick,monotone,monthDays,movStats,mtitle,multEventChart,multLines,na.delete,na.detail.response,na.include,na.keep,na.pattern,na.retain,naclus,nafitted.delete,Names2names,naplot,napredict.delete,naprint.delete,naprint.keep,naresid.delete,naresid.keep,nCoincident,nFm,nmChoice,nobsY,nomiss,nstr,num.denom.setup,num.intercepts,numeric.string,numericScale,oPar,optionsCmds,ordGridFun,ordTestpo,outerText,pairUpDiff,panel.bpplot,panel.Dotplot,panel.Ecdf,panel.plsmo,panel.xYplot,parGrid,partition.matrix,partition.vector,pasteFit,pBlock,pc1,plotCorrM,plotCorrPrecision,plotp,plotlyM,plotlyParm,plotlySave,plotmathTranslate,plotMultSim,plotpsummaryM,plsmo,pMedian,pngNeedle,pomodm,popower,posamsize,prepanel.Dotplot,prepanel.Ecdf,prepanel.xYplot,princmp,printL,print.char.matrix,prList,prn,propsPO,propsTrans,prselect,prType,pstamp,putHcap,putHfig,putKey,putKeyEmpty,qcrypt,Quantile,Quantile2,r2describe,R2Measures,rcorr,rcorr.cens,rcorrcens,rcorrp.cens,rcspline.eval,rcspline.plot,rcspline.restate,rcsplineFunction,read.xportDataload,readSAScsv,redun,reformM,reLabelled,rendHTML,replace.substring.wild,reShape,responseSummary,restoreHattrib,rlegend,rlegendg,rm.boot,rMultinom,roundN,roundPOSIXt,runifChanged,runParallel,samplesize.bin,sas.codes,sas.get,sas.get.macro,sasdsLabels,sasxport.get,Save,scat1d,score.binary,sedit,sepUnitsTrans,seqFreq,setParNro,setTrellis,show.col,show.pch,showPsfrag,simMarkovOrd,simplifyDims,simPOcuts,simRegOrd,sKey,smean.cl.boot,smean.cl.normal,smean.sd,smean.sdl,smearingEst,smedian.hilow,solvet,somers2,soprobMarkovOrd,soprobMarkovOrdm,spearman,spearman.test,spearman2,spikecomp,spower,spss.get,src,stat_plsmo,stata.get,StatPlsmo,stepfun.eval,stratify,strgraphwrap,string.bounding.box,string.break.line,stringDims,stripChart,subplot,substi,substi.source,substring.location,substring2,summarize,summaryD,summaryDp,summaryM,summaryP,summaryRc,summaryS,symbol.freq,sys,t.test.cluster,table_formatpct,table_freq,table_latexdefs,table_N,table_pc,table_trio,tabulr,termsDrop,testCharDateTime,testDateTime,tex,tobase64image,transace,transcan,translate,trap.rule,trellis.strip.blank,truncPOSIXt,uncbind,units,unPaste,upData,upFirst,valueLabel,valueName,valueTags,valueUnit,var.inner,varclus,vlab,Weibull2,whichClosek,whichClosePW,whichClosest,wtd.Ecdf,wtd.loess.noiter,wtd.mean,wtd.quantile,wtd.rank,wtd.table,wtd.var,xInch,xless,xy.group,xYplot,xySortNoDupNoNA,yearDays,yInch,ynbind,zoom,"[<-.discrete","consolidate<-","is.na<-.discrete","label<-","length<-.discrete","substring2<-","valueLabel<-","valueName<-","valueTags<-","valueUnit<-") useDynLib(Hmisc, .registration=TRUE, .fixes="F_") diff --git a/NEWS b/NEWS index 1aea794..9f540dc 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,10 @@ +Changes in version 5.2-0 (2024-10-19) + * added requireNamespace for safer in qcrypt + * require R 4.2.0 or later, to handle pipes + * new function pMedian for pseudomedian (Hodges-Lehmann one-sample estimator) + * added pMedian to describe function + * Converted all Ratfor code to Fortran 2018 + Changes in version 5.1-3 (2024-05-18) * areg.boot: intercepted use of xlim when y is categorical; need to improve logic on xout= * redun: stored scores matrix in returned object, suitable for r2describe @@ -7,7 +14,7 @@ Changes in version 5.1-3 (2024-05-18) * qcrypt: added service=NA to bypass keyring * runParallel: changed progress directory explicitly to /tmp for Mac * printL: new function - * redun: got around model.frame problem so that negated variables are omitted from analysis + * redun: got around model.frame problem so that negated variables are omitted from analysis; also fixed printing of missing values Changes in version 5.1-2 (2024-03-09) * dotchartpl: remove length(ugroup) criterion (unsure what this was doing and made lower and upper be not plotted) diff --git a/R/describe.s b/R/describe.s index 7d3daa1..04dc8db 100644 --- a/R/describe.s +++ b/R/describe.s @@ -153,10 +153,11 @@ describe.vector <- function(x, descript, exclude.missing=TRUE, digits=4, } else counts <- c(counts, format(sum(weights * x) / sum(weights), ...)) lab <- c(lab, "Mean") - if(! weighted) { - gmd <- format(GiniMd(xnum), ...) - counts <- c(counts, gmd) - lab <- c(lab, "Gmd") + if(! weighted & n.unique > 2) { + pmedian <- format(pMedian(xnum)) + gmd <- format(GiniMd(xnum), ...) + counts <- c(counts, pmedian, gmd) + lab <- c(lab, "pMedian", "Gmd") } } else if(n.unique == 1) { counts <- c(counts, format(x.unique)) @@ -1384,6 +1385,7 @@ html_describe_con <- function(x, sparkwidth=200, Distinct = a('distinct'), Info = a('Info'), Mean = k['Mean'], + pMedian = k['pMedian'], Gmd = k['Gmd'], ' ' = '' , # place for sparkline '.05' = r('.05'), @@ -1443,7 +1445,7 @@ html_describe_con <- function(x, sparkwidth=200, gt::tab_header(title=gt::md(title), subtitle=subtitle) |> gt::tab_style(style=gt::cell_text(align='center'), locations=gt::cells_column_labels( - columns=c(n, Missing, Distinct, Info, Mean, Gmd))) |> + columns=c(n, Missing, Distinct, Info, Mean, pMedian, Gmd))) |> gt::tab_style(style=gt::cell_text(size='small'), locations=gt::cells_body(columns=Label)) |> gt::text_transform(locations=gt::cells_body(columns=' '), @@ -1539,6 +1541,7 @@ html_describe_cat <- function(x, w=200, freq=c('chart', 'table'), Info = a('Info'), Sum = a('Sum'), Mean = p('Mean'), + pMedian = p('pMedian'), Gmd = p('Gmd'), tab = tab, type = type, @@ -1553,11 +1556,11 @@ html_describe_cat <- function(x, w=200, freq=c('chart', 'table'), spik <- a$tab[a$type == 'spark'] a$tab[a$type == 'spark'] <- '' - for(i in c('Units', 'Format', 'Sum', 'Mean', 'Info', 'Gmd')) + for(i in c('Units', 'Format', 'Sum', 'Mean', 'Info', 'pMedian', 'Gmd')) if(all(is.na(a[[i]])) || all(a[[i]] == '')) a[[i]] <- NULL center_cols <- intersect(names(a), - c('n', 'Missing', 'Distinct', 'Info', 'Sum', 'Mean', 'Gmd')) + c('n', 'Missing', 'Distinct', 'Info', 'Sum', 'Mean', 'pMedian', 'Gmd')) b <- gt::gt(a) |> gt::tab_header(title=gt::md(title), subtitle=subtitle) |> @@ -1577,13 +1580,14 @@ html_describe_cat <- function(x, w=200, freq=c('chart', 'table'), gt::cols_label(tab = ' ') |> gt::sub_missing(missing_text='') + if('Units' %in% names(a)) b <- b |> gt::tab_style(style=gt::cell_text(size='small', font='arial'), locations=gt::cells_body(columns=Units)) if('Gmd' %in% names(a)) b <- b |> gt::cols_label(Gmd ~ gt::html("Gini\u2009|\u394|")) - +saveRDS(b, '/tmp/b.rds') b } diff --git a/R/pMedian.r b/R/pMedian.r new file mode 100644 index 0000000..2da0195 --- /dev/null +++ b/R/pMedian.r @@ -0,0 +1,27 @@ +#' Pseudomedian +#' +#' Uses fast Fortran code to compute the pseudomedian of a numeric vector. The pseudomedian is the median of all possible midpoints of two observations. The pseudomedian is also called the Hodges-Lehmann one-sample estimator. The Fortran code is was originally from JF Monahan, and was converted to C++ in the `DescTools` package. It has been converted to Fortran 2018 here. +#' @title pMedian +#' @param x a numeric vector +#' @param na.rm set to `TRUE` to exclude `NA`s before computing the pseudomedian +#' +#' @return a scalar numeric value +#' @export +#' @md +#' @seealso , +#' @examples +#' x <- c(1:4, 10000) +#' pMedian(x) +#' # Compare with brute force calculation and with wilcox.test +#' w <- outer(x, x, '+') +#' median(w[lower.tri(w, diag=TRUE)]) / 2 +#' wilcox.test(x, conf.int=TRUE) + +pMedian <- function(x, na.rm = FALSE) { + if(na.rm) x <- x[! is.na(x)] + n <- length(x) + if(n == 0) return(NA_real_) + if(n == 1) return(as.double(x)) + .Fortran(F_hlqest, as.double(sort(x)), as.double(runif(1000)), as.integer(n), result=double(1))$result + +} diff --git a/R/qcrypt.r b/R/qcrypt.r index a9f22da..f86d73d 100644 --- a/R/qcrypt.r +++ b/R/qcrypt.r @@ -48,6 +48,9 @@ qcrypt <- function(obj, base, service='R-keyring-service', file) { stop('you must install the getPass package to use qcrypt') if(! requireNamespace('qs', quietly=TRUE)) stop('you must install the qs package to use qcrypt') +if(! requireNamespace('safer', quietly=TRUE)) + stop('you must install the safer package to use qcrypt') + if(is.na(service)) { prompt <- if(! missing(base)) 'Define password for storing encrypted data: ' diff --git a/inst/tests/pMedian.r b/inst/tests/pMedian.r new file mode 100644 index 0000000..90757b3 --- /dev/null +++ b/inst/tests/pMedian.r @@ -0,0 +1,94 @@ +require(Hmisc) + +pm <- function(x, method=c('uniroot', 'optimize', 'bruteforce', 'cpp', 'fortran')) { + method <- match.arg(method) + if(method == 'fortran') return(pMedian(x)) + if(method == 'cpp') return(hlqest(x)) + if(method == 'bruteforce') { # good for N < 1000 but still slower than otheres + w <- outer(x, x, '+') + return(median(w[lower.tri(w, diag=TRUE)]) / 2) + } + r <- range(x) + # Compute [0.1] normalized ranks just for numeric reasons + nrank <- function(x) (rank(x) - 1) / (length(x) - 1) + if(method == 'uniroot') { + obj <- function(mu) sum(nrank(abs(x - mu)) * sign(x - mu)) + xs <- seq(r[1], r[2], length=100) + o <- sapply(xs, obj) + zero <- approx(o, xs, xout=0)$y + R <- uniroot(obj, r, tol=1e-6) + if(abs(obj(zero)) < abs(R$f.root)) return(zero) + R$root + } else { + obj <- function(mu) sum(nrank(abs(x - mu)) * sign(x - mu)) ^ 2 + optimize(obj, r, tol=1e-7)$minimum + # optim(r[2], obj, method=if(FALSE)'L-BFGS-B' else 'Brent', lower=r[1], upper=r[3])$par + } +} + +# Compile C++ code from DescTools package +hl <- readLines('https://raw.githubusercontent.com/AndriSignorell/DescTools/refs/heads/master/src/hlqest.cpp') +Rcpp::sourceCpp(code=paste(hl, collapse='\n')) +# Code is also stored in ~/rr/pseudomedian/hlqest.cpp + +set.seed(1) +w <- NULL +for(n in seq(5, 1000, by=5)) { + x <- round(runif(n) * 100) + w <- rbind(w, data.frame(n=n, + brute=pm(x, 'bruteforce'), cpp=pm(x, 'cpp'), fortran=pm(x, 'fortran'), + root=pm(x, 'uniroot'), mean=mean(x), median=median(x))) +} +for(i in 2: 4) for(j in (i + 1) : 5) cat(i, j, max(abs(w[,i] - w[,j])), '\n') +with(w, plot(n, root - cpp)) +table(w$cpp >= pmin(w$mean, w$median) & w$cpp <= pmax(w$mean, w$median)) +w[which.max(abs(w$cpp - w$root)),] +set.seed(1) +x <- round(runif(5) * 100) +x +x <- c(20, 27, 37, 57, 91) +xcpp <- pm(x, 'cpp') +xroot <- pm(x, 'uniroot') +pm(x, 'optimize') +c(xcpp, xroot) +xs <- seq(20, 91, by=0.1) +obj <- sapply(xs, function(z) sum(rank(abs(x - z)) * sign(x - z))) +approx(obj, xs, xout=0) +plot(xs, obj) +abline(v=xcpp, col='red') +abline(v=xroot, col='blue') +abline(h=0, col=gray(0.5)) + +g <- function(method, n=10000) { + for(i in 1 : 1000) { + x <- runif(n) + pm(x, method) + } +} + +system.time(g('bruteforce', 200)) +system.time(g('uniroot', 1000)) +system.time(g('cpp')) +system.time(g('fortran')) + +x <- runif(5e6) +system.time(a <- pm(x, 'cpp')) +system.time(b <- pm(x, 'fortran')) +a; b + +set.seed(1) +N <- integer(1000); dif <- numeric(1000) +for(i in 1 : 1000) { + n <- sample(1 : 75000, 1) + x <- round(runif(n) * 100) + N[i] <- n + dif[i] <- abs(pm(x, 'cpp') - pm(x, 'fortran')) +} +hist(dif, nclass=50) +plot(N, dif) + + +# Original Fortran function fails at N=46507 +# The square root of the largest integer in R is (2^31 - 1) is 46341 +# So the problem is in the computation of N*(N+1)/2 pairs in hlqest +# C++ function uses long integers; changed Fortran to use INTEGER*8 likewise diff --git a/man/Hmisc-internal.Rd b/man/Hmisc-internal.Rd index ca31e71..bb4674e 100644 --- a/man/Hmisc-internal.Rd +++ b/man/Hmisc-internal.Rd @@ -83,6 +83,7 @@ \alias{F_rcorr} \alias{F_wclosepw} \alias{F_wclosest} +\alias{F_hlqest} \description{Internal Hmisc functions.} \details{These are not to be called by the user or are undocumented.} diff --git a/man/describe.Rd b/man/describe.Rd index 9904375..31d2d31 100644 --- a/man/describe.Rd +++ b/man/describe.Rd @@ -115,8 +115,8 @@ Note: As discussed in Cox and Longton (2008), Stata Technical Bulletin 8(4) pp. 557, the term "unique" has been replaced with "distinct" in the output (but not in parameter names). -When \code{weights} are not used, Gini's mean difference is computed for -numeric variables. This is a robust measure of dispersion that is the +When \code{weights} are not used, the pseudomedian and Gini's mean difference are computed for +numeric variables. The pseudomedian is labeled \code{pMedian} and is the median of all possible pairwise averages. It is a robust and efficient measure of location that equals the mean and median for symmetric distributions. It is also called the Hodges-Lehmann one-sample estimator. Gini's mean difference is a robust measure of dispersion that is the mean absolute difference between any pairs of observations. In simple output Gini's difference is labeled \code{Gmd}. @@ -378,7 +378,7 @@ Vanderbilt University } \seealso{ \code{\link{spikecomp}}, \code{\link{sas.get}}, \code{\link{quantile}}, -\code{\link{GiniMd}}, +\code{\link{GiniMd}}, \code{\link{pMedian}}, \code{\link{table}}, \code{\link{summary}}, \code{\link{model.frame.default}}, \code{\link{naprint}}, \code{\link{lapply}}, \code{\link{tapply}}, diff --git a/man/latex.Rd b/man/latex.Rd index 1296488..7e7040b 100644 --- a/man/latex.Rd +++ b/man/latex.Rd @@ -40,7 +40,7 @@ commands that are similar to those in the \code{S.to.latex} procedure in the \code{s.to.latex} package (Chambers and Hastie, 1993). \code{latex.function} can also produce \code{verbatim} output or output that works with the \code{Sweavel} -LaTeX style at \url{https://biostat.app.vumc.org/wiki/Main/SweaveTemplate}. +LaTeX style. \code{latex.list} calls \code{latex} recursively for each element in the argument. @@ -462,7 +462,7 @@ dvigv(object, \dots) \code{center="none"} to use no centering. \code{centerline} can be useful when objects besides a \code{tabular} are enclosed in a single \code{table} environment. - This option was implemented by Markus Jäntti + This option was implemented by Markus J�ntti \email{markus.jantti@iki.fi} of Abo Akademi University. } \item{landscape}{ diff --git a/man/pMedian.Rd b/man/pMedian.Rd new file mode 100644 index 0000000..d176305 --- /dev/null +++ b/man/pMedian.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pMedian.r +\name{pMedian} +\alias{pMedian} +\title{pMedian} +\usage{ +pMedian(x, na.rm = FALSE) +} +\arguments{ +\item{x}{a numeric vector} + +\item{na.rm}{set to \code{TRUE} to exclude \code{NA}s before computing the pseudomedian} +} +\value{ +a scalar numeric value +} +\description{ +Pseudomedian +} +\details{ +Uses fast Fortran code to compute the pseudomedian of a numeric vector. The pseudomedian is the median of all possible midpoints of two observations. The pseudomedian is also called the Hodges-Lehmann one-sample estimator. The Fortran code is was originally from JF Monahan, and was converted to C++ in the \code{DescTools} package. It has been converted to Fortran 2018 here. +} +\examples{ +x <- c(1:4, 10000) +pMedian(x) +# Compare with brute force calculation and with wilcox.test +w <- outer(x, x, '+') +median(w[lower.tri(w, diag=TRUE)]) / 2 +wilcox.test(x, conf.int=TRUE) +} +\seealso{ +\url{https://dl.acm.org/doi/10.1145/1271.319414}, \url{https://www4.stat.ncsu.edu/~monahan/jul10/} +} diff --git a/man/popower.Rd b/man/popower.Rd index ddbca98..926643a 100644 --- a/man/popower.Rd +++ b/man/popower.Rd @@ -28,19 +28,19 @@ differences in mean or median on the original scale. a PO model to demonstrate the natural sampling variability that causes estimated odds ratios to vary over cutoffs of Y. -\code{propsPO} uses \code{\link{ggplot2}} to plot a stacked bar chart of +\code{propsPO} uses \code{ggplot2} to plot a stacked bar chart of proportions stratified by a grouping variable (and optionally a stratification variable), with an optional additional graph showing what the proportions would be had proportional odds held and an odds ratio was applied to the proportions in a reference group. If the result is passed to \code{ggplotly}, customized tooltip hover text will appear. -\code{propsTrans} uses \code{\link{ggplot2}} to plot all successive +\code{propsTrans} uses \code{ggplot2} to plot all successive transition proportions. \code{formula} has the state variable on the left hand side, the first right-hand variable is time, and the second right-hand variable is a subject ID variable.\ -\code{multEventChart} uses \code{\link{ggplot2}} to plot event charts +\code{multEventChart} uses \code{ggplot2} to plot event charts showing state transitions, account for absorbing states/events. It is based on code written by Lucy D'Agostino McGowan posted at \url{https://livefreeordichotomize.com/posts/2020-05-21-survival-model-detective-1/}. diff --git a/man/transcan.Rd b/man/transcan.Rd index 6c2f575..632e2b6 100644 --- a/man/transcan.Rd +++ b/man/transcan.Rd @@ -360,7 +360,7 @@ invertTabulated(x, y, freq=rep(1,length(x)), slower approach is to use \code{impcat="multinom"} to fit a multinomial logistic model to the categorical variable, at the last iteraction of the - \code{transcan} algorithm. This uses the \code{\link{multinom}} + \code{transcan} algorithm. This uses the \code{\link[nnet]{multinom}} function in the \pkg{nnet} library of the \pkg{MASS} package (which is assumed to have been installed by the user) to fit a polytomous logistic model to the current working transformations of all the diff --git a/src/Hmisc.h b/src/Hmisc.h index 0dd80dd..df8c187 100644 --- a/src/Hmisc.h +++ b/src/Hmisc.h @@ -7,6 +7,14 @@ #include #include "R_ext/Error.h" +#ifdef __cplusplus +#define EXTERNC extern "C" +#else +#define EXTERNC +#endif + +#undef EXTERNC + #ifdef _SPLUS_ # define STRING_ELT(x,i) (CHARACTER_POINTER(x)[i]) # define TO_CHAR(x) (x) diff --git a/src/cidxcn.f b/src/cidxcn.f deleted file mode 100644 index 0b121c4..0000000 --- a/src/cidxcn.f +++ /dev/null @@ -1,88 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 -C------------------------------------------------------------------------ -C Compute c-index (c) and Brown-Hollander-Krowar-Goodman-Kruskal-Somer -C rank correlation (gamma) between X and Y with censoring indicator E -C Also returns number of relevant, concordant, and uncertain pairs -C (nrel, nconc, nuncert) and estimated s.d. of gamma (sd) using -C Quade formula (see SAS PROC MATPAR). Pairs with tied x are -C excluded if outx=.TRUE. -C -C F. Harrell 27Nov90 -C Modification of SAS Procedure KGKC (1980) -C------------------------------------------------------------------------- - subroutine cidxcn(x,y,e,n,nrel,nconc,nuncert,c,gamma,sd,outx) - implicit double precision (a-h,o-z) - double precision x(n),y(n),dx,dy - logical e(n),outx - double precision nrel,nuncert,nconc - nconc=0d0 - nrel=0d0 - nuncert=0d0 - sumr=0d0 - sumr2=0d0 - sumw=0d0 - sumw2=0d0 - sumrw=0d0 - do23000 i=1,n - wi=0d0 - ri=0d0 - do23002 j=1,n - if(j.ne.i)then - dx=x(i)-x(j) - dy=y(i)-y(j) - if(dx.ne.0. .or. .not.outx)then - if((e(i).and.dy.lt.0.).or.(e(i).and..not.e(j).and.dy.eq.0.))then - if(dx.lt.0.)then - nconc=nconc+1d0 - wi=wi+1d0 - else - if(dx.eq.0.)then - nconc=nconc+.5d0 - else - wi=wi-1d0 - endif - endif - nrel=nrel+1d0 - ri=ri+1d0 - else - if((e(j).and.dy.gt.0.).or.(e(j).and..not.e(i).and.dy.eq.0.))then - if(dx.gt.0.)then - nconc=nconc+1d0 - wi=wi+1d0 - else - if(dx.eq.0.)then - nconc=nconc+.5d0 - else - wi=wi-1d0 - endif - endif - nrel=nrel+1d0 - ri=ri+1d0 - else - if(.not.(e(i).and.e(j)))then - nuncert=nuncert+1d0 - endif - endif - endif - endif - endif -23002 continue -23003 continue - sumr=sumr+ri - sumr2=sumr2+ri*ri - sumw=sumw+wi - sumw2=sumw2+wi*wi - sumrw=sumrw+ri*wi -23000 continue -23001 continue - c=nconc/nrel - gamma=2.*(c-.5) -Ccall dblepr('sumr',4,sumr,1) -Ccall dblepr('sumw',4,sumw,1) -Ccall dblepr('sumr2',5,sumr2,1) -Ccall dblepr('sumw2',5,sumw2,1) -Ccall dblepr('sumrw',5,sumrw,1) - sd=sumr2*sumw**2-2d0*sumr*sumw*sumrw+sumw2*sumr**2 - sd=2.*dsqrt(sd)/sumr/sumr - return - end diff --git a/src/cidxcn.f90 b/src/cidxcn.f90 new file mode 100644 index 0000000..13d5c5c --- /dev/null +++ b/src/cidxcn.f90 @@ -0,0 +1,95 @@ +! Converted from Ratfor to Fortran 2018 by Claude AI 2024-10-26 +subroutine cidxcn(x, y, e, n, nrel, nconc, nuncert, c, gamma, sd, outx) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + + ! Arguments + integer, intent(in) :: n + real(real64), intent(in) :: x(n), y(n) + logical, intent(in) :: e(n), outx + real(real64), intent(out) :: nrel, nconc, nuncert, c, gamma, sd + + ! Local variables + integer :: i, j + real(real64) :: dx, dy, wi, ri + real(real64) :: sumr, sumr2, sumw, sumw2, sumrw + + ! Initialize counters and sums + nconc = 0.0_real64 + nrel = 0.0_real64 + nuncert = 0.0_real64 + sumr = 0.0_real64 + sumr2 = 0.0_real64 + sumw = 0.0_real64 + sumw2 = 0.0_real64 + sumrw = 0.0_real64 + + ! Main computation loop + do i = 1, n + wi = 0.0_real64 + ri = 0.0_real64 + + do j = 1, n + if (j /= i) then + dx = x(i) - x(j) + dy = y(i) - y(j) + + ! Check if pair should be considered based on x ties + if (dx /= 0.0_real64 .or. (.not. outx)) then + ! Case 1: i is uncensored and y(i) < y(j) + ! or i is uncensored, j is censored, and y values are equal + if ((e(i) .and. dy < 0.0_real64) .or. & + (e(i) .and. .not. e(j) .and. dy == 0.0_real64)) then + + if (dx < 0.0_real64) then + nconc = nconc + 1.0_real64 + wi = wi + 1.0_real64 + else if (dx == 0.0_real64) then + nconc = nconc + 0.5_real64 + else + wi = wi - 1.0_real64 + end if + nrel = nrel + 1.0_real64 + ri = ri + 1.0_real64 + + ! Case 2: j is uncensored and y(j) > y(i) + ! or j is uncensored, i is censored, and y values are equal + else if ((e(j) .and. dy > 0.0_real64) .or. & + (e(j) .and. .not. e(i) .and. dy == 0.0_real64)) then + + if (dx > 0.0_real64) then + nconc = nconc + 1.0_real64 + wi = wi + 1.0_real64 + else if (dx == 0.0_real64) then + nconc = nconc + 0.5_real64 + else + wi = wi - 1.0_real64 + end if + nrel = nrel + 1.0_real64 + ri = ri + 1.0_real64 + + ! Case 3: Uncertain pair (both censored) + else if (.not. (e(i) .and. e(j))) then + nuncert = nuncert + 1.0_real64 + end if + end if + end if + end do + + ! Accumulate sums for standard deviation calculation + sumr = sumr + ri + sumr2 = sumr2 + ri * ri + sumw = sumw + wi + sumw2 = sumw2 + wi * wi + sumrw = sumrw + ri * wi + end do + + ! Calculate c-index and gamma + c = nconc / nrel + gamma = 2.0_real64 * (c - 0.5_real64) + + ! Calculate standard deviation using Quade formula + sd = sumr2 * sumw**2 - 2.0_real64 * sumr * sumw * sumrw + sumw2 * sumr**2 + sd = 2.0_real64 * sqrt(sd) / (sumr * sumr) + +end subroutine cidxcn diff --git a/src/cidxcp.f b/src/cidxcp.f deleted file mode 100644 index be8bf7c..0000000 --- a/src/cidxcp.f +++ /dev/null @@ -1,92 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine cidxcp(x1,x2,y,e,n,method,outx, nrel,nuncer,c1,c2,gamma - *1,gamma2,gamma,sd,c12,c21) - implicit double precision (a-h,o-z) - double precision x1(n),x2(n),y(n) - logical e(n),outx - integer n,method,i,j - double precision nrel,nuncer,nconc1,nconc2,c1,c2,gamma1,gamma2, su - *mr,sumr2,sumw,sumw2,sumrw, wi,ri,sumc,c12,c21,gamma,sd - double precision dx,dx2,dy - nconc1=0d0 - nconc2=0d0 - nrel=0d0 - nuncer=0d0 - sumr=0d0 - sumr2=0d0 - sumw=0d0 - sumw2=0d0 - sumrw=0d0 - sumc=0d0 - do23000 i=1,n - wi=0d0 - ri=0d0 - do23002 j=1,n - dx=x1(i)-x1(j) - dx2=x2(i)-x2(j) - if((i.ne.j) .and. (.not.outx .or. dx.ne.0. .or. dx2.ne.0.))then - dy=y(i)-y(j) - if((e(i).and.(dy.lt.0.)).or.(e(i).and..not.e(j).and.(dy.eq.0.)))th - *en - nrel=nrel+1d0 - nconc1=nconc1+(z(dx.lt.0.)+.5d0*z(dx.eq.0.)) - nconc2=nconc2+(z(dx2.lt.0.)+.5d0*z(dx2.eq.0.)) - ri=ri+1d0 - if(method.eq.1)then - wi=wi+(z(dx.lt.dx2)-z(dx.gt.dx2)) - sumc=sumc+z(dx.lt.dx2) - else - wi=wi+(z(dx.lt.0..and.dx2.ge.0.)-z(dx.gt.0..and.dx2.le.0.)) - sumc=sumc+z(dx.lt.0..and.dx2.ge.0.) - endif - else - if((e(j).and.(dy.gt.0.)).or.(e(j).and..not.e(i).and.(dy.eq.0.)))th - *en - nrel=nrel+1d0 - nconc1=nconc1+(z(dx.gt.0.)+.5d0*z(dx.eq.0.)) - nconc2=nconc2+(z(dx2.gt.0.)+.5d0*z(dx2.eq.0.)) - ri=ri+1d0 - if(method.eq.1)then - wi=wi+(z(dx.gt.dx2)-z(dx.lt.dx2)) - sumc=sumc+z(dx.gt.dx2) - else - wi=wi+(z(dx.gt.0..and.dx2.le.0.)-z(dx.lt.0..and.dx2.ge.0.)) - sumc=sumc+z(dx.gt.0..and.dx2.le.0.) - endif - else - if(.not.(e(i).and.e(j)))then - nuncer=nuncer+1d0 - endif - endif - endif - endif -23002 continue -23003 continue - sumr=sumr+ri - sumr2=sumr2+ri*ri - sumw=sumw+wi - sumw2=sumw2+wi*wi - sumrw=sumrw+ri*wi -23000 continue -23001 continue - c1=nconc1/nrel - gamma1=2d0*(c1-.5d0) - c2=nconc2/nrel - gamma2=2d0*(c2-.5d0) - gamma=sumw/sumr - sd=sumr2*sumw**2-2d0*sumr*sumw*sumrw+sumw2*sumr**2 - sd=2d0*dsqrt(sd)/sumr/sumr - c12=sumc/sumr - c21=sumc/sumr-gamma - return - end - function z(a) - double precision z - logical a - if(a)then - z=1d0 - else - z=0d0 - endif - return - end diff --git a/src/cidxcp.f90 b/src/cidxcp.f90 new file mode 100644 index 0000000..9b0c453 --- /dev/null +++ b/src/cidxcp.f90 @@ -0,0 +1,94 @@ +! Converted from Ratfor to Fortran 2018 by ChatGPT 2024-10-26 +! Instruction to ChatGPT: +! Convert the following Ratfor code to Fortran 2018 using Fortran ISO environment without using module, and make the z function strictly internal to the subroutine + +subroutine cidxcp(x1, x2, y, e, n, method, outx, nrel, nuncer, c1, c2, gamma1, gamma2, gamma, sd, c12, c21) + use ISO_FORTRAN_ENV, only: REAL64, INT32 + implicit none + integer(INT32), intent(in) :: n, method + real(REAL64), intent(in) :: x1(n), x2(n), y(n) + logical, intent(in) :: e(n), outx + real(REAL64), intent(out) :: nrel, nuncer, c1, c2, gamma1, gamma2, gamma, sd, c12, c21 + + integer(INT32) :: i, j + real(REAL64) :: nconc1, nconc2, sumr, sumr2, sumw, sumw2, sumrw, sumc + real(REAL64) :: wi, ri, dx, dx2, dy + + ! Initialize accumulation variables + nconc1 = 0.0_REAL64; nconc2 = 0.0_REAL64; nrel = 0.0_REAL64; nuncer = 0.0_REAL64 + sumr = 0.0_REAL64; sumr2 = 0.0_REAL64; sumw = 0.0_REAL64 + sumw2 = 0.0_REAL64; sumrw = 0.0_REAL64; sumc = 0.0_REAL64 + + do i = 1, n + wi = 0.0_REAL64; ri = 0.0_REAL64 + do j = 1, n + dx = x1(i) - x1(j) + dx2 = x2(i) - x2(j) + + if ((i /= j) .and. (.not. outx .or. dx /= 0.0_REAL64 .or. dx2 /= 0.0_REAL64)) then + dy = y(i) - y(j) + + if ((e(i) .and. (dy < 0.0_REAL64)) .or. (e(i) .and. .not. e(j) .and. (dy == 0.0_REAL64))) then + nrel = nrel + 1.0_REAL64 + nconc1 = nconc1 + (z(dx < 0.0_REAL64) + 0.5_REAL64 * z(dx == 0.0_REAL64)) + nconc2 = nconc2 + (z(dx2 < 0.0_REAL64) + 0.5_REAL64 * z(dx2 == 0.0_REAL64)) + ri = ri + 1.0_REAL64 + + if (method == 1) then + wi = wi + (z(dx < dx2) - z(dx > dx2)) + sumc = sumc + z(dx < dx2) + else + wi = wi + (z(dx < 0.0_REAL64 .and. dx2 >= 0.0_REAL64) - z(dx > 0.0_REAL64 .and. dx2 <= 0.0_REAL64)) + sumc = sumc + z(dx < 0.0_REAL64 .and. dx2 >= 0.0_REAL64) + endif + + else if ((e(j) .and. (dy > 0.0_REAL64)) .or. (e(j) .and. .not. e(i) .and. (dy == 0.0_REAL64))) then + nrel = nrel + 1.0_REAL64 + nconc1 = nconc1 + (z(dx > 0.0_REAL64) + 0.5_REAL64 * z(dx == 0.0_REAL64)) + nconc2 = nconc2 + (z(dx2 > 0.0_REAL64) + 0.5_REAL64 * z(dx2 == 0.0_REAL64)) + ri = ri + 1.0_REAL64 + + if (method == 1) then + wi = wi + (z(dx > dx2) - z(dx < dx2)) + sumc = sumc + z(dx > dx2) + else + wi = wi + (z(dx > 0.0_REAL64 .and. dx2 <= 0.0_REAL64) - z(dx < 0.0_REAL64 .and. dx2 >= 0.0_REAL64)) + sumc = sumc + z(dx > 0.0_REAL64 .and. dx2 <= 0.0_REAL64) + endif + + else if (.not. (e(i) .and. e(j))) then + nuncer = nuncer + 1.0_REAL64 + endif + endif + end do + + sumr = sumr + ri + sumr2 = sumr2 + ri * ri + sumw = sumw + wi + sumw2 = sumw2 + wi * wi + sumrw = sumrw + ri * wi + end do + + ! Calculate outputs + c1 = nconc1 / nrel + gamma1 = 2.0_REAL64 * (c1 - 0.5_REAL64) + c2 = nconc2 / nrel + gamma2 = 2.0_REAL64 * (c2 - 0.5_REAL64) + gamma = sumw / sumr + sd = 2.0_REAL64 * sqrt(sumr2 * sumw**2 - 2.0_REAL64 * sumr * sumw * sumrw + sumw2 * sumr**2) / sumr / sumr + c12 = sumc / sumr + c21 = sumc / sumr - gamma + + contains + + real(REAL64) function z(a) + logical, intent(in) :: a + if (a) then + z = 1.0_REAL64 + else + z = 0.0_REAL64 + endif + end function z + + end subroutine cidxcp + \ No newline at end of file diff --git a/src/hlqest.f90 b/src/hlqest.f90 new file mode 100644 index 0000000..dac4c9c --- /dev/null +++ b/src/hlqest.f90 @@ -0,0 +1,144 @@ + subroutine hlqest(x, ran, nn, retval) + use ISO_FORTRAN_ENV + implicit none + integer(int32), intent(in) :: nn + real(real64), intent(in) :: x(nn) + real(real64), intent(in) :: ran(1000) + real(real64), intent(out) :: retval + + integer(int64) :: n, np, sm, l, k1, k2, sq, mdlrow, k, i, j, iran + real (real64) :: am, amn, amx + integer(int64), allocatable :: lb(:), rb(:), q(:) + + n = nn + iran = 0 + + ! Handle small n + if (n < 3) then + retval = (x(1) + x(n)) / 2.0 + return + endif + + ! Allocate memory for lb, rb, and q + allocate(lb(n), rb(n), q(n)) + + ! Total number of pairs and the median(s) index + np = n * (n + 1) / 2 + k1 = (np + 1) / 2 + k2 = (np + 2) / 2 + + ! Initialize lb and rb + lb = [(i, i=1, n)] ! Vectorized initialization of lb + rb = n ! Vectorized initialization of rb to n for all elements + + sm = np ! Total number of pairs + l = 0 ! Number of pairs less than those in set at step m + am = x((n + 1) / 2) + x((n + 2) / 2) + + ! Main iterative process + do + j = n + + ! Vectorized loop to compute q and sq + q = 0 + sq = 0 + + do i = 1, n + do while (j >= i) + if (x(i) + x(j) < am) then + q(i) = j - i + 1 + sq = sq + q(i) + exit + endif + j = j - 1 + end do + end do + + ! Handle case where partitions are the same (ties) + if (sq == l) then + amx = x(1) + x(1) + amn = x(n) + x(n) + do i = 1, n + if (lb(i) > rb(i)) cycle + amn = min(amn, x(lb(i)) + x(i)) ! Vectorized min within bounds + amx = max(amx, x(rb(i)) + x(i)) ! Vectorized max within bounds + end do + am = (amn + amx) / 2.0 + if ((am <= amn) .or. (am > amx)) am = amx + if ((amn /= amx) .and. (sm /= 2)) cycle + retval = am / 2.0 + exit + endif ! end if(sq == l) + + ! Handle if we're nearly done (values on the border) + if (sq == k2 - 1 .or. sq == k1) then + amx = x(1) + x(1) + amn = x(n) + x(n) + + ! Loop over all rows to find max of those < am and min of those >= am, considering q + do i = 1, n + if (q(i) > 0) amx = max(amx, x(i) + x(i + q(i) - 1)) + ! There are q(i) pairs in row i that satisfy x[i] + x[j] < am + if (q(i) < (n - i + 1)) amn = min(amn, x(i) + x(i + q(i))) + ! There are (n - i + 1 - q(i)) pairs in row i that satisfy x[i] + x[j] >= am + end do + + ! Handle case where k1 < k2 (when we need the average of two values) + if (k1 < k2) then + retval = (amn + amx) / 4.0 + else if (sq == k1) then + retval = amx / 2.0 + else if (sq == k1 - 1) then + retval = amn / 2.0 + else + retval = (amn + amx) / 4.0 + end if + exit + endif ! end if(sq == k2 - 1 .or ...) + + ! Update bounds based on the size of sq using vectorization + if (sq < k1) then + lb = [(i + q(i), i=1, n)] ! Vectorized update of left bounds + else + rb = [(i + q(i) - 1, i=1, n)] ! Vectorized update of right bounds + endif + + ! Recalculate sm and l using vectorized array operations + l = sum(lb - [(i, i=1, n)]) ! Vectorized summation of lb - i + sm = sum(rb - lb + 1) ! Vectorized summation of pairs still in set + + ! If we still have more than 2 pairs, pick a random row and continue + if (sm > 2) then + iran = iran + 1 + if(iran > 1000) iran = 1 + k = int(ran(iran) * dble(sm)) + do i = 1, n + j = i + if (k <= rb(i) - lb(i)) exit + k = k - rb(i) + lb(i) - 1 + end do + mdlrow = (lb(j) + rb(j)) / 2 + am = x(j) + x(mdlrow) + cycle + endif + + amx = x(1) + x(1) ! ? + amn = x(n) + x(n) ! ? + + ! Otherwise, finalize and exit + do i = 1, n + if (lb(i) > rb(i)) cycle + amn = min(amn, x(lb(i)) + x(i)) + amx = max(amx, x(rb(i)) + x(i)) + end do + am = (amx + amn) / 2.0 + if ((am <= amn) .or. (am > amx)) am = amx + if ((am /= amx .and. sm /= 2)) cycle + retval = am / 2.0 + exit + end do + + ! Deallocate memory + deallocate(lb, rb, q) + return + end subroutine hlqest diff --git a/src/hoeffd.f b/src/hoeffd.f deleted file mode 100644 index ef41eac..0000000 --- a/src/hoeffd.f +++ /dev/null @@ -1,120 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine hoeffd(xx, n, p, dmat, aadmat, madmat, npair, x, y, rx, - * ry, rj) - implicit double precision (a-h,o-z) - integer p, npair(p,p) - double precision xx(n,p), dmat(p,p), aadmat(p,p), madmat(p,p), x(n - *), y(n), rx(n), ry(n), rj(n), maxad - do23000 i=1, p - np=0 - do23002 k=1, n - if(xx(k,i) .lt. 1d49)then - np = np + 1 - endif -23002 continue -23003 continue - npair(i,i) = np - do23006 j=(i+1),p - m = 0 - do23008 k=1,n - xk = xx(k,i) - yk = xx(k,j) - if(xk .lt. 1d49 .and. yk .lt. 1d49)then - m = m + 1 - x(m) = xk - y(m) = yk - endif -23008 continue -23009 continue - npair(i,j) = m - if(m .gt. 4)then - call hoeff(x, y, m, d, aad, maxad, rx, ry, rj) - dmat(i,j) = d - aadmat(i,j) = aad - madmat(i,j) = maxad - else - dmat(i,j) = 1d50 - endif -23006 continue -23007 continue -23000 continue -23001 continue - do23014 i=1,p - dmat(i,i) = 1d0/30d0 - do23016 j=(i+1),p - dmat(j,i) = dmat(i,j) - npair(j,i) = npair(i,j) - aadmat(j,i) = aadmat(i,j) - madmat(j,i) = madmat(i,j) -23016 continue -23017 continue -23014 continue -23015 continue - return - end - subroutine hoeff(x, y, n, d, aad, maxad, rx, ry, rj) - implicit double precision (a-h,o-z) - double precision x(n), y(n), rx(n), ry(n), rj(n), maxad - call jrank(x, y, n, rx, ry, rj) - q = 0d0 - r = 0d0 - s = 0d0 - aad = 0d0 - maxad = 0d0 - z = n - do23018 i=1,n - rxi = rx(i) - ryi = ry(i) - rji = rj(i) - ad = dabs((rji/z) - (rxi/z)*(ryi/z)) - aad = aad + ad - maxad = dmax1(maxad, ad) - q = q + (rxi-1d0)*(rxi-2d0)*(ryi-1d0)*(ryi-2d0) - r = r + (rxi-2d0)*(ryi-2d0)*(rji-1d0) - s = s + (rji-1d0)*(rji-2d0) -23018 continue -23019 continue - aad = aad / z - d = (q-2d0*(z-2d0)*r+(z-2d0)*(z-3d0)*s)/z/(z-1d0)/(z-2d0)/(z-3d0)/ - *(z-4d0) - return - end - subroutine jrank(x, y, n, rx, ry, r) - integer n - double precision x(n), y(n), rx(n), ry(n), r(n), cx, cy, ri, rix, - *riy, xi, yi - do23020 i=1,n - xi = x(i) - yi = y(i) - ri = 1d0 - rix = 1d0 - riy = 1d0 - do23022 j=1,n - if(i .ne. j)then - cx = 0d0 - if(x(j) .lt. xi)then - cx = 1d0 - endif - if(x(j) .eq. xi)then - cx = .5d0 - endif - cy = 0d0 - if(y(j) .lt. yi)then - cy = 1d0 - endif - if(y(j) .eq. yi)then - cy = .5d0 - endif - rix = rix + cx - riy = riy + cy - ri = ri + cx*cy - endif -23022 continue -23023 continue - rx(i) = rix - ry(i) = riy - r(i) = ri -23020 continue -23021 continue - return - end diff --git a/src/hoeffd.f90 b/src/hoeffd.f90 new file mode 100644 index 0000000..7e64753 --- /dev/null +++ b/src/hoeffd.f90 @@ -0,0 +1,146 @@ +! Converted from Ratfor to Fortran 2018 by Claude AI 2024-10-26 +! +! Hoeffding's D Statistics Implementation +! Modern Fortran without module structure + +subroutine hoeffd(xx, n, p, dmat, aadmat, madmat, npair, x, y, rx, ry, rj) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + + ! Arguments + integer, intent(in) :: n, p + real(real64), intent(in) :: xx(n,p) + real(real64), intent(out) :: dmat(p,p), aadmat(p,p), madmat(p,p) + integer, intent(out) :: npair(p,p) + real(real64), intent(out) :: x(n), y(n), rx(n), ry(n), rj(n) + + ! Local variables + integer :: i, j, k, m, np + real(real64) :: xk, yk, d, aad, maxad + real(real64), parameter :: LARGE_VAL = 1.0e49_real64 + + ! Initialize arrays + dmat = 0.0_real64 + aadmat = 0.0_real64 + madmat = 0.0_real64 + npair = 0 + + do i = 1, p + np = count(xx(:,i) < LARGE_VAL) + npair(i,i) = np + + do j = i+1, p + m = 0 + do k = 1, n + xk = xx(k,i) + yk = xx(k,j) + if (xk < LARGE_VAL .and. yk < LARGE_VAL) then + m = m + 1 + x(m) = xk + y(m) = yk + end if + end do + + npair(i,j) = m + if (m > 4) then + call hoeff(x, y, m, d, aad, maxad, rx, ry, rj) + dmat(i,j) = d + aadmat(i,j) = aad + madmat(i,j) = maxad + else + dmat(i,j) = LARGE_VAL + end if + end do + end do + + ! Fill lower triangle + do i = 1, p + dmat(i,i) = 1.0_real64/30.0_real64 + do j = i+1, p + dmat(j,i) = dmat(i,j) + npair(j,i) = npair(i,j) + aadmat(j,i) = aadmat(i,j) + madmat(j,i) = madmat(i,j) + end do + end do + +end subroutine hoeffd + +subroutine hoeff(x, y, n, d, aad, maxad, rx, ry, rj) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + + ! Arguments + integer, intent(in) :: n + real(real64), intent(in) :: x(n), y(n) + real(real64), intent(out) :: d, aad, maxad + real(real64), intent(out) :: rx(n), ry(n), rj(n) + + ! Local variables + integer :: i + real(real64) :: q, r, s, z, rxi, ryi, rji, ad + + call jrank(x, y, n, rx, ry, rj) + + q = 0.0_real64 + r = 0.0_real64 + s = 0.0_real64 + aad = 0.0_real64 + maxad = 0.0_real64 + z = real(n, real64) + + do i = 1, n + rxi = rx(i) + ryi = ry(i) + rji = rj(i) + ad = abs((rji/z) - (rxi/z)*(ryi/z)) + aad = aad + ad + maxad = max(maxad, ad) + q = q + (rxi-1.0_real64)*(rxi-2.0_real64)*(ryi-1.0_real64)*(ryi-2.0_real64) + r = r + (rxi-2.0_real64)*(ryi-2.0_real64)*(rji-1.0_real64) + s = s + (rji-1.0_real64)*(rji-2.0_real64) + end do + + aad = aad / z + d = (q-2.0_real64*(z-2.0_real64)*r+(z-2.0_real64)*(z-3.0_real64)*s) / & + (z*(z-1.0_real64)*(z-2.0_real64)*(z-3.0_real64)*(z-4.0_real64)) + +end subroutine hoeff + +subroutine jrank(x, y, n, rx, ry, r) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + + ! Arguments + integer, intent(in) :: n + real(real64), intent(in) :: x(n), y(n) + real(real64), intent(out) :: rx(n), ry(n), r(n) + + ! Local variables + integer :: i, j + real(real64) :: xi, yi, ri, rix, riy, cx, cy + + do i = 1, n + xi = x(i) + yi = y(i) + ri = 1.0_real64 + rix = 1.0_real64 + riy = 1.0_real64 + + do j = 1, n + if (i /= j) then + cx = merge(0.5_real64, merge(1.0_real64, 0.0_real64, x(j) < xi), x(j) == xi) + cy = merge(0.5_real64, merge(1.0_real64, 0.0_real64, y(j) < yi), y(j) == yi) + + rix = rix + cx + riy = riy + cy + ri = ri + cx*cy + end if + end do + + rx(i) = rix + ry(i) = riy + r(i) = ri + end do + +end subroutine jrank diff --git a/src/init.c b/src/init.c index 944972b..fb9c133 100644 --- a/src/init.c +++ b/src/init.c @@ -22,6 +22,7 @@ extern void F77_NAME(rcorr)(void *, void *, void *, void *, void *, void *, void /* extern void F77_NAME(wcidxy)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); */ extern void F77_NAME(wclosepw)(void *, void *, void *, void *, void *, void *, void *, void *); extern void F77_NAME(wclosest)(void *, void *, void *, void *, void *); +extern void F77_NAME(hlqest)(void *, void *, void *, void *); static const R_CallMethodDef CallEntries[] = { {"do_mchoice_match", (DL_FUNC) &do_mchoice_match, 3}, @@ -40,6 +41,7 @@ static const R_FortranMethodDef FortranEntries[] = { /* {"wcidxy", (DL_FUNC) &F77_NAME(wcidxy), 11}, */ {"wclosepw", (DL_FUNC) &F77_NAME(wclosepw), 8}, {"wclosest", (DL_FUNC) &F77_NAME(wclosest), 5}, + {"hlqest", (DL_FUNC) &F77_NAME(hlqest), 4}, {NULL, NULL, 0} }; diff --git a/src/jacklins.f b/src/jacklins.f deleted file mode 100644 index f33d710..0000000 --- a/src/jacklins.f +++ /dev/null @@ -1,23 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine jacklins(x, w, n, k, res) - integer n, k, l - double precision x(n), w(n-1,k), res(n,k), sj - do23000 l = 1, k - do23002 j = 1, n - sj = 0d0 - do23004 i = 1, n - if(i .lt. j)then - sj = sj + w(i,l) * x(i) - endif - if(i .gt. j)then - sj = sj + w(i-1,l) * x(i) - endif -23004 continue -23005 continue - res(j,l) = sj -23002 continue -23003 continue -23000 continue -23001 continue - return - end diff --git a/src/jacklins.f90 b/src/jacklins.f90 new file mode 100644 index 0000000..ad5bdf6 --- /dev/null +++ b/src/jacklins.f90 @@ -0,0 +1,31 @@ +! Converted from Ratfor 2024-10-26 using ChatGPT and the following request: +! Convert the following Ratfor code to Fortran 2018 using Fortran ISO environment without using module + +subroutine jacklins(x, w, n, k, res) + use ISO_FORTRAN_ENV, only: REAL64, INT32 + implicit none + + integer(INT32), intent(in) :: n, k + real(REAL64), intent(in) :: x(n), w(n-1, k) + real(REAL64), intent(out) :: res(n, k) + integer(INT32) :: l, j, i + real(REAL64) :: sj + + ! Compute leave-out-one linear statistics for each column + do l = 1, k + do j = 1, n + sj = 0.0_REAL64 + do i = 1, n + if (i < j) then + sj = sj + w(i, l) * x(i) + else if (i > j) then + sj = sj + w(i - 1, l) * x(i) + endif + end do + res(j, l) = sj + end do + end do + + return + end subroutine jacklins + diff --git a/src/maxempr.f b/src/maxempr.f deleted file mode 100644 index af9a5da..0000000 --- a/src/maxempr.f +++ /dev/null @@ -1,68 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine maxempr(ax, ay, x, y, n, w, h, z, area, rect) - implicit double precision (a-h,o-z) - integer n - double precision ax(2), ay(2), x(n), y(n), z(3), rect(4), maxr, li - maxr = z(1) * dabs(ay(2) - ay(1)) - rect(1) = z(2) - rect(2) = ay(1) - rect(3) = z(3) - rect(4) = ay(2) - do23000 i=1,n - tl = ax(1) - tr = ax(2) - if(i .lt. n)then - do23004 j=(i+1),n - if(x(j) .gt. tl .and. x(j) .lt. tr)then - area = (tr - tl) * (y(j) - y(i)) - if(area .gt. maxr .and. ((tr - tl) .gt. w) .and. ((y(j) - y(i)) .g - *t. h))then - maxr = area - rect(1) = tl - rect(2) = y(i) - rect(3) = tr - rect(4) = y(j) - endif - if(x(j) .gt. x(i))then - tr = x(j) - else - tl = x(j) - endif - endif -23004 continue -23005 continue - endif - area = (tr - tl) * (ay(2) - y(i)) - if(area .gt. maxr .and. ((tr - tl) .gt. w) .and. ((ay(2) - y(i)) . - *gt. h))then - maxr = area - rect(1) = tl - rect(2) = y(i) - rect(3) = tr - rect(4) = ay(2) - endif - ri = ax(2) - li = ax(1) - do23014 k=1,n - if(y(k) .lt. y(i) .and. x(k) .gt. x(i))then - ri = dmin1(ri, x(k)) - endif - if(y(k) .lt. y(i) .and. x(k) .lt. x(i))then - li = dmax1(li, x(k)) - endif -23014 continue -23015 continue - area = (ri - li) * (ay(2) - y(i)) - if(area .gt. maxr .and. ((ri - li) .gt. w) .and. ((y(i) - ay(1)) . - *gt. h))then - maxr = area - rect(1) = li - rect(2) = ay(1) - rect(3) = ri - rect(4) = y(i) - endif -23000 continue -23001 continue - area = maxr - return - end diff --git a/src/maxempr.f90 b/src/maxempr.f90 new file mode 100644 index 0000000..3b560ca --- /dev/null +++ b/src/maxempr.f90 @@ -0,0 +1,81 @@ +! Converted from Ratfor 2024-10-26 using ChatGPT and the following request: +! Convert the following Ratfor code to Fortran 2018 using Fortran ISO environment without using module + +subroutine maxempr(ax, ay, x, y, n, w, h, z, area, rect) + use ISO_FORTRAN_ENV, only: REAL64, INT32 + implicit none + + integer(INT32), intent(in) :: n + real(REAL64), intent(in) :: ax(2), ay(2), x(n), y(n), w, h, z(3) + real(REAL64), intent(out) :: area, rect(4) + real(REAL64) :: maxr, tl, tr, li, ri + real(REAL64) :: area_candidate + integer(INT32) :: i, j, k + + ! Initialize max area and set initial rectangle boundaries + maxr = z(1) * abs(ay(2) - ay(1)) + rect(1) = z(2) + rect(2) = ay(1) + rect(3) = z(3) + rect(4) = ay(2) + + ! Iterate through points to determine largest rectangles + do i = 1, n + tl = ax(1) + tr = ax(2) + + if (i < n) then + do j = i + 1, n + if (x(j) > tl .and. x(j) < tr) then + ! Check horizontal slices (j == i + 1) + area_candidate = (tr - tl) * (y(j) - y(i)) + if (area_candidate > maxr .and. (tr - tl) > w .and. (y(j) - y(i)) > h) then + maxr = area_candidate + rect(1) = tl + rect(2) = y(i) + rect(3) = tr + rect(4) = y(j) + end if + + ! Update boundaries based on current point + if (x(j) > x(i)) then + tr = x(j) + else + tl = x(j) + end if + end if + end do + end if + + ! Check open rectangles above (x(i), y(i)) + area_candidate = (tr - tl) * (ay(2) - y(i)) + if (area_candidate > maxr .and. (tr - tl) > w .and. (ay(2) - y(i)) > h) then + maxr = area_candidate + rect(1) = tl + rect(2) = y(i) + rect(3) = tr + rect(4) = ay(2) + end if + + ! Check open rectangles below (x(i), y(i)) + ri = ax(2) + li = ax(1) + do k = 1, n + if (y(k) < y(i) .and. x(k) > x(i)) ri = min(ri, x(k)) + if (y(k) < y(i) .and. x(k) < x(i)) li = max(li, x(k)) + end do + + area_candidate = (ri - li) * (y(i) - ay(1)) + if (area_candidate > maxr .and. (ri - li) > w .and. (y(i) - ay(1)) > h) then + maxr = area_candidate + rect(1) = li + rect(2) = ay(1) + rect(3) = ri + rect(4) = y(i) + end if + end do + + area = maxr + return + end subroutine maxempr + \ No newline at end of file diff --git a/src/rcorr.f b/src/rcorr.f deleted file mode 100644 index 16c5302..0000000 --- a/src/rcorr.f +++ /dev/null @@ -1,106 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine rcorr(xx, n, p, itype, dmat, npair, x, y, rx, ry, work, - * iwork) - integer p, npair(p,p) - double precision xx(n,p), dmat(p,p), x(n), y(n), rx(n), ry(n), wor - *k(n) - integer iwork(n) - double precision sumx, sumx2, sumy, sumy2, sumxy, z, a, b, xk, yk, - * d - sumx = 0d0 - sumy = 0d0 - sumx2 = 0d0 - sumy2 = 0d0 - sumxy = 0d0 - do23000 i = 1, p - np = 0 - do23002 k = 1, n - if(xx(k, i) .lt. 1d49)then - np = np + 1 - endif -23002 continue -23003 continue - npair(i, i) = np - do23006 j = (i + 1), p - m = 0 - if(itype .eq. 1)then - sumx = 0d0 - sumy = 0d0 - sumx2 = 0d0 - sumy2 = 0d0 - sumxy = 0d0 - endif - do23010 k = 1, n - xk = xx(k, i) - yk = xx(k, j) - if(xk .lt. 1d49 .and. yk .lt. 1d49)then - m = m + 1 - if(itype .eq. 1)then - a = xk - b = yk - sumx = sumx + a - sumx2 = sumx2 + a * a - sumy = sumy + b - sumy2 = sumy2 + b * b - sumxy = sumxy + a * b - else - x(m) = xk - y(m) = yk - endif - endif -23010 continue -23011 continue - npair(i, j) = m - if(m .gt. 1)then - if(itype .eq. 1)then - z = m - d = (sumxy - sumx * sumy / z) / dsqrt((sumx2 - sumx * sumx / z) * - *(sumy2 - sumy * sumy / z)) - else - call docorr(x, y, m, d, rx, ry, work, iwork) - endif - dmat(i, j) = d - else - dmat(i, j) = 1d50 - endif -23006 continue -23007 continue -23000 continue -23001 continue - do23020 i = 1, p - dmat(i, i) = 1d0 - do23022 j = (i + 1), p - dmat(j, i) = dmat(i, j) - npair(j, i) = npair(i, j) -23022 continue -23023 continue -23020 continue -23021 continue - return - end - subroutine docorr(x, y, n, d, rx, ry, work, iwork) - double precision x(1), y(1), rx(n), ry(n), work(1) - integer iwork(1) - double precision sumx, sumx2, sumy, sumy2, sumxy, a, b, z, d - call rank(n, x, work, iwork, rx) - call rank(n, y, work, iwork, ry) - sumx = 0d0 - sumx2 = 0d0 - sumy = 0d0 - sumy2 = 0d0 - sumxy = 0d0 - do23024 i = 1, n - a = rx(i) - b = ry(i) - sumx = sumx + a - sumx2 = sumx2 + a * a - sumy = sumy + b - sumy2 = sumy2 + b * b - sumxy = sumxy + a * b -23024 continue -23025 continue - z = n - d = (sumxy - sumx * sumy / z) / dsqrt((sumx2 - sumx * sumx / z) * - *(sumy2 - sumy * sumy / z)) - return - end diff --git a/src/rcorr.f90 b/src/rcorr.f90 new file mode 100644 index 0000000..5785d31 --- /dev/null +++ b/src/rcorr.f90 @@ -0,0 +1,123 @@ +! Converted from Ratfor 2024-10-26 using ChatGPT and the following request: +! Convert the following Ratfor code to Fortran 2018 using Fortran ISO environment without using module + +subroutine rcorr(xx, n, p, itype, dmat, npair, x, y, rx, ry, work, iwork) + use ISO_FORTRAN_ENV + implicit none + integer, intent(in) :: n, p, itype + real(REAL64), intent(in) :: xx(n, p) + real(REAL64), intent(out) :: dmat(p, p) + integer, intent(out) :: npair(p, p) + real(REAL64), intent(out) :: x(n), y(n), rx(n), ry(n), work(n) + integer, intent(out) :: iwork(n) + real(REAL64) :: sumx, sumx2, sumy, sumy2, sumxy, z, a, b, xk, yk, d + integer :: i, j, k, m, np + + ! Initialize sums to prevent warnings + sumx = 0.0_real64 + sumy = 0.0_real64 + sumx2 = 0.0_real64 + sumy2 = 0.0_real64 + sumxy = 0.0_real64 + + ! Loop through each pair of variables + do i = 1, p + np = 0 + do k = 1, n + if (xx(k, i) < 1.0e49_real64) np = np + 1 + end do + npair(i, i) = np + + do j = i + 1, p + m = 0 + if (itype == 1) then + sumx = 0.0_real64 + sumy = 0.0_real64 + sumx2 = 0.0_real64 + sumy2 = 0.0_real64 + sumxy = 0.0_real64 + end if + + do k = 1, n + xk = xx(k, i) + yk = xx(k, j) + if (xk < 1.0e49_real64 .and. yk < 1.0e49_real64) then + m = m + 1 + if (itype == 1) then + a = xk + b = yk + sumx = sumx + a + sumx2 = sumx2 + a * a + sumy = sumy + b + sumy2 = sumy2 + b * b + sumxy = sumxy + a * b + else + x(m) = xk + y(m) = yk + end if + end if + end do + + npair(i, j) = m + if (m > 1) then + if (itype == 1) then + z = real(m, REAL64) + d = (sumxy - sumx * sumy / z) / sqrt((sumx2 - sumx * sumx / z) * (sumy2 - sumy * sumy / z)) + else + call docorr(x, y, m, d, rx, ry, work, iwork) + end if + dmat(i, j) = d + else + dmat(i, j) = 1.0e50_real64 + end if + end do + end do + + ! Set the diagonal and symmetry + do i = 1, p + dmat(i, i) = 1.0_real64 + do j = i + 1, p + dmat(j, i) = dmat(i, j) + npair(j, i) = npair(i, j) + end do + end do + + return + end subroutine rcorr + + subroutine docorr(x, y, n, d, rx, ry, work, iwork) + use ISO_FORTRAN_ENV + implicit none + integer, intent(in) :: n + real(REAL64), intent(in) :: x(n), y(n) + real(REAL64), intent(out) :: d + real(REAL64), intent(out) :: rx(n), ry(n), work(n) + integer, intent(out) :: iwork(n) + real(REAL64) :: sumx, sumx2, sumy, sumy2, sumxy, a, b, z + integer :: i + + call rank(n, x, work, iwork, rx) + call rank(n, y, work, iwork, ry) + + sumx = 0.0_real64 + sumx2 = 0.0_real64 + sumy = 0.0_real64 + sumy2 = 0.0_real64 + sumxy = 0.0_real64 + + do i = 1, n + a = rx(i) + b = ry(i) + sumx = sumx + a + sumx2 = sumx2 + a * a + sumy = sumy + b + sumy2 = sumy2 + b * b + sumxy = sumxy + a * b + end do + + z = real(n, REAL64) + d = (sumxy - sumx * sumy / z) / sqrt((sumx2 - sumx * sumx / z) * (sumy2 - sumy * sumy / z)) + + return + end subroutine docorr + \ No newline at end of file diff --git a/src/wclosest.f b/src/wclosest.f deleted file mode 100644 index 1c85c02..0000000 --- a/src/wclosest.f +++ /dev/null @@ -1,57 +0,0 @@ -C Output from Public domain Ratfor, version 1.01 - subroutine wclosest(w, x, lw, lx, j) - implicit double precision (a-h,o-z) - integer lw, lx, j(lw) - double precision w(lw), x(lx) - do23000 i=1,lw - wi=w(i) - dmin=1d40 - m=0 - do23002 k=1,lx - d = dabs(x(k) - wi) - if(d .lt. dmin)then - dmin = d - m = k - endif -23002 continue -23003 continue - j(i) = m -23000 continue -23001 continue - return - end - subroutine wclosepw(w, x, r, f, lw, lx, xd, j) - implicit double precision (a-h,o-z) - double precision w(lw),x(lx),r(lw),xd(lx) - integer lw, lx, j(lw) - do23006 i=1, lw - wi = w(i) - dmean = 0d0 - do23008 k=1, lx - xd(k) = dabs(x(k) - wi) - dmean = dmean + xd(k) -23008 continue -23009 continue - dmean = f * dmean / lx - sump = 0d0 - do23010 k=1, lx - z = min(xd(k)/dmean, 1d0) - xd(k) = (1d0 - z**3)**3 - sump = sump + xd(k) -23010 continue -23011 continue - prob = 0d0 - ri = r(i) - m = 1 - do23012 k=1, lx - prob = prob + xd(k) / sump - if(ri .gt. prob)then - m = m + 1 - endif -23012 continue -23013 continue - j(i) = m -23006 continue -23007 continue - return - end diff --git a/src/wclosest.f90 b/src/wclosest.f90 new file mode 100644 index 0000000..8bcd885 --- /dev/null +++ b/src/wclosest.f90 @@ -0,0 +1,72 @@ +! Converted from Ratfor 2024-10-26 using ChatGPT and the following request: +! Convert the following Ratfor code to Fortran 2018 using Fortran ISO environment without using module + +subroutine wclosest(w, x, lw, lx, j) + use ISO_FORTRAN_ENV, only: REAL64 + implicit none + integer, intent(in) :: lw, lx + integer, intent(out) :: j(lw) + real(REAL64), intent(in) :: w(lw), x(lx) + real(REAL64) :: wi, dmin, d + integer :: i, k, m + + do i = 1, lw + wi = w(i) + dmin = 1.0e40_real64 + m = 0 + do k = 1, lx + d = abs(x(k) - wi) + if (d < dmin) then + dmin = d + m = k + end if + end do + j(i) = m + end do + + return + end subroutine wclosest + + subroutine wclosepw(w, x, r, f, lw, lx, xd, j) + use ISO_FORTRAN_ENV, only: REAL64 + implicit none + integer, intent(in) :: lw, lx + integer, intent(out) :: j(lw) + real(REAL64), intent(in) :: w(lw), x(lx), r(lw), f + real(REAL64), intent(out) :: xd(lx) + real(REAL64) :: wi, dmean, sump, z, prob, ri + integer :: i, k, m + + do i = 1, lw + wi = w(i) + dmean = 0.0_real64 + + ! Calculate mean distance + do k = 1, lx + xd(k) = abs(x(k) - wi) + dmean = dmean + xd(k) + end do + dmean = f * dmean / real(lx, REAL64) + + ! Adjust probabilities + sump = 0.0_real64 + do k = 1, lx + z = min(xd(k) / dmean, 1.0_real64) + xd(k) = (1.0_real64 - z**3)**3 + sump = sump + xd(k) + end do + + ! Determine index based on probabilities + prob = 0.0_real64 + ri = r(i) + m = 1 + do k = 1, lx + prob = prob + xd(k) / sump + if (ri > prob) m = m + 1 + end do + j(i) = m + end do + + return + end subroutine wclosepw + \ No newline at end of file