[r-cran-raschsampler] 03/19: Import Upstream version 0.8-4
Andreas Tille
tille at debian.org
Mon Dec 12 10:45:17 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-raschsampler.
commit e6628aa55651a41adb82550fc55069350ad0a3f6
Author: Andreas Tille <tille at debian.org>
Date: Mon Dec 12 11:20:43 2016 +0100
Import Upstream version 0.8-4
---
COPYING | 340 +++++++++++++++++++++++++++
COPYRIGHTS | 14 ++
DESCRIPTION | 13 ++
NAMESPACE | 10 +
R/phi.range.R | 9 +
R/rsampler.R | 48 ++++
R/rsctrl.R | 15 ++
R/rserror.R | 32 +++
R/rsextrmat.R | 15 ++
R/rsextrobj.R | 32 +++
R/rstats.R | 25 ++
R/rsunpack.R | 35 +++
R/summary.RSctr.R | 13 ++
R/summary.RSmpl.R | 16 ++
R/summary.RSmplext.R | 15 ++
data/xmpl.rda | Bin 0 -> 3239 bytes
data/xmplbig.rda | Bin 0 -> 131982 bytes
man/RSctr.Rd | 45 ++++
man/RSmpl.Rd | 46 ++++
man/RaschSampler.package.Rd | 80 +++++++
man/phi.range.Rd | 25 ++
man/rsampler.Rd | 63 +++++
man/rsctrl.Rd | 86 +++++++
man/rsextrmat.Rd | 28 +++
man/rsextrobj.Rd | 67 ++++++
man/rstats.Rd | 66 ++++++
man/summary.RSctr.Rd | 19 ++
man/summary.RSmpl.Rd | 31 +++
man/xmpl.Rd | 25 ++
src/RaschSampler.f90 | 550 ++++++++++++++++++++++++++++++++++++++++++++
30 files changed, 1763 insertions(+)
diff --git a/COPYING b/COPYING
new file mode 100755
index 0000000..fbdd65f
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,340 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/COPYRIGHTS b/COPYRIGHTS
new file mode 100755
index 0000000..a962cbe
--- /dev/null
+++ b/COPYRIGHTS
@@ -0,0 +1,14 @@
+COPYRIGHT STATUS for RaschSampler
+---------------------------------
+
+This R code is
+
+ Copyright (C) 2006 Reinhold Hatzinger, Patrick Mair and Norman Verhelst
+
+The Fortran code is
+
+ Copyright (C) 2006 Norman Verhelst
+
+All code is subject to the GNU General Public License, Version 2. See
+the file COPYING for the exact conditions under which you may
+redistribute it.
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100755
index 0000000..4fe204a
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,13 @@
+Package: RaschSampler
+Type: Package
+Title: Rasch Sampler
+Version: 0.8-4
+Date: 2010-03-22
+Author: Reinhold Hatzinger, Patrick Mair, Norman Verhelst
+Maintainer: <reinhold.hatzinger at wu.ac.at>
+Depends: stats
+Description: Sampling binary matrices with fixed margins
+License: GPL
+Repository: CRAN
+Date/Publication: 2010-03-23 07:40:32
+Packaged: 2010-03-22 17:33:01 UTC; hatz
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100755
index 0000000..8ba3487
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,10 @@
+useDynLib(RaschSampler)
+export(rsampler)
+export(rstats)
+export(rsctrl)
+export(rsextrobj)
+export(rsextrmat)
+export(phi.range)
+S3method(summary, RSctr)
+S3method(summary, RSmpl)
+S3method(summary, RSmplext)
diff --git a/R/phi.range.R b/R/phi.range.R
new file mode 100755
index 0000000..fe0dbbe
--- /dev/null
+++ b/R/phi.range.R
@@ -0,0 +1,9 @@
+"phi.range" <-
+function(mat){
+ cmat<-cor(mat)+diag(NA,ncol(mat))
+ ma<-max(cmat,na.rm=TRUE)
+ mi<-min(cmat,na.rm=TRUE)
+ RET <- ma-mi
+ RET
+}
+
diff --git a/R/rsampler.R b/R/rsampler.R
new file mode 100755
index 0000000..05e30a3
--- /dev/null
+++ b/R/rsampler.R
@@ -0,0 +1,48 @@
+"rsampler" <-
+function(inpmat,controls=rsctrl()){
+
+ if (!(class(controls)=="RSctr"))
+ stop("controls is not a control object - see help(\"rsctrl\")")
+
+ n <- dim(inpmat)[1]
+ k <- dim(inpmat)[2]
+ burn_in <- controls$burn_in
+ n_eff <- controls$n_eff
+ step <- controls$step
+ seed <- controls$seed
+ tfixed <- controls$tfixed
+
+ if (seed == 0) {
+ # generates random seed in the range [536870911,772830910]
+ seed <- as.integer(as.double(format(Sys.time(), "%H%M%OS3"))*1000)
+ + 2**29 - 1
+ }
+
+ # allocation of memory for simulated matrices
+ vec<-vector( length = (n_eff+1)*n*trunc((k+31)/32) )
+ ier<-0
+
+ # calls the external Fortran subroutine sampler
+ # simulated matrices are returned in vec
+ RET<-.Fortran("sampler",
+ n=as.integer(n),
+ k=as.integer(k),
+ inpmat=as.integer(inpmat),
+ tfixed=as.logical(tfixed),
+ burn_in=as.integer(burn_in),
+ n_eff=as.integer(n_eff),
+ step=as.integer(step),
+ seed=as.integer(seed),
+ outvec=as.integer(vec),
+ ier=as.integer(ier)
+ )
+ n_tot <- n_eff+1
+ if (RET$ier>0) {
+ rserror(RET$ier)
+ } else {
+ RET<-c(RET[1:8],n_tot=n_eff+1,RET[9:10])
+ class(RET)<-"RSmpl"
+ RET
+ }
+}
+
diff --git a/R/rsctrl.R b/R/rsctrl.R
new file mode 100755
index 0000000..57d10ca
--- /dev/null
+++ b/R/rsctrl.R
@@ -0,0 +1,15 @@
+"rsctrl" <-
+function(burn_in=100, n_eff=100, step=16,
+ seed=0, tfixed=FALSE)
+{ ier <- 0
+ if(n_eff < 0 | n_eff > 8191) ier <- ier + 4
+ if(burn_in < 0) ier <- ier + 8
+ if(step <= 0) ier <- ier + 16
+ if(seed < 0 | seed > (2**31-2)) ier <- ier + 128
+ if (ier>0) rserror(ier)
+ RET<-list(burn_in=burn_in, n_eff=n_eff, step=step,
+ seed=seed, tfixed=tfixed)
+ class(RET)<-"RSctr"
+ RET
+}
+
diff --git a/R/rserror.R b/R/rserror.R
new file mode 100755
index 0000000..0bec7e5
--- /dev/null
+++ b/R/rserror.R
@@ -0,0 +1,32 @@
+"rserror" <-
+function(err)
+{
+ bin2int<-function(x){
+ r<-vector(mode="integer",length=0)
+ while(x>1){
+ r<-c(x%%2,r)
+ x<-floor(x/2)
+ }
+ r<-c(1,r)
+ rev(r)
+ }
+ errortxt<-c("\tn < 0 or n > 4096\n",
+ "\tk < 0 or k > 128\n",
+ "\tn_eff < 0 or n_eff > 8191\n",
+ "\tburn_in < 0\n",
+ "\tstep <= 0\n",
+ "\tone or more entries in the input matrix are different from 1 or 0\n",
+ "\tthe input matrix is of Guttman form; the sample space has only one element\n",
+ "\tseed < 0 or seed > 2147483646\n"
+ )
+
+ if (err == 0) {
+ cat("\nno error\n")
+ } else {
+ x <- bin2int(err)
+ errstring <- paste("\n",paste(errortxt[(1:length(x))*x],sep="",collapse=""))
+ stop(errstring, call.=FALSE)
+ }
+ invisible(err)
+}
+
diff --git a/R/rsextrmat.R b/R/rsextrmat.R
new file mode 100755
index 0000000..1e80397
--- /dev/null
+++ b/R/rsextrmat.R
@@ -0,0 +1,15 @@
+"rsextrmat" <-
+function(RSobj, mat.no = 1)
+{
+ obj.name <- deparse(substitute(RSobj))
+ if (!(class(RSobj)=="RSmpl" || class(RSobj)=="RSmplext")){
+ err.text<-paste(obj.name," is not a sample object - see help(\"rsextrobj\")",sep ="",collapse="")
+ stop(err.text)
+ }
+ if(mat.no > RSobj$n_tot)
+ stop("\n\tElement ",mat.no," not available (",obj.name," has ", RSobj$n_tot, " elements).")
+ obj<-rsextrobj(RSobj, start = mat.no, end = mat.no)
+ RET<-rstats(obj, function(x) matrix(x, nr = obj$n))[[1]]
+ RET
+}
+
diff --git a/R/rsextrobj.R b/R/rsextrobj.R
new file mode 100755
index 0000000..0081713
--- /dev/null
+++ b/R/rsextrobj.R
@@ -0,0 +1,32 @@
+"rsextrobj" <-
+function(RSobj,start=1,end=8192)
+{
+ obj.name <- deparse(substitute(RSobj))
+ if (!(class(RSobj)=="RSmpl" || class(RSobj)=="RSmplext")){
+ err.text<-paste(obj.name,"not a sample object - see help(\"rsextrobj\")",sep ="",collapse="")
+ stop(err.text)
+ }
+
+ n_tot <- RSobj$n_tot
+ if (end>n_tot) end<-n_tot
+ n <- RSobj$n
+ k <- RSobj$k
+ nwords <- c(trunc((k+31)/32))
+
+ objnew <- RSobj
+ l_one_mat <- n*nwords
+ b <- (start-1)*l_one_mat+1
+ e <- end*l_one_mat
+ objnew$outvec <- RSobj$outvec[b:e]
+ objnew$n_tot <- end-start+1
+ if (start==1) {
+ objnew$n_eff <- objnew$n_tot - 1
+ } else {
+ objnew$n_eff <- objnew$n_tot
+ }
+ class(objnew)="RSmplext"
+
+ RET<-objnew
+ RET
+}
+
diff --git a/R/rstats.R b/R/rstats.R
new file mode 100755
index 0000000..b771856
--- /dev/null
+++ b/R/rstats.R
@@ -0,0 +1,25 @@
+"rstats" <-
+function(RSobj,userfunc,...)
+{
+ obj.name <- deparse(substitute(RSobj))
+ if (!(class(RSobj)=="RSmpl" || class(RSobj)=="RSmplext")){
+ err.text<-paste(obj.name," is not a sample object - see help(\"rsextrobj\")",sep ="",collapse="")
+ stop(err.text)
+ }
+
+ # extracts simulated matrices into three dimensional array sim
+ n_tot <- RSobj$n_tot
+ n <- RSobj$n
+ k <- RSobj$k
+ nwords <- c(trunc((k+31)/32))
+
+ # store coded simulated matrices in list with n_eff+1 elements
+ sim<-split(RSobj$outvec,gl(n_tot,n*nwords))
+
+
+ # decode simulated matrices and apply user function
+ #RET<-unlist(lapply(sim,rsunpack,n,k,nwords,userfunc))
+ RET<-lapply(sim,rsunpack,n,k,nwords,userfunc,...)
+ RET
+}
+
diff --git a/R/rsunpack.R b/R/rsunpack.R
new file mode 100755
index 0000000..32fe668
--- /dev/null
+++ b/R/rsunpack.R
@@ -0,0 +1,35 @@
+"rsunpack" <-
+function(x,n,k,nwords,userfunc,...){
+ # check for NAs (-2^31 not defined in R) in simulated matrices
+ # set value to 0
+
+ nas<-FALSE
+ if (k>=32) {
+ idx<-(1:length(x))[is.na(x)] # indexvector for NAs
+ nas<-(length(idx)>0)
+ x[idx]<-0
+ }
+
+ t<-vector(length=n*k)
+
+ # calls unpacking routine
+ out<-.Fortran("unpack",
+ as.integer(x),
+ as.integer(nwords),
+ mat=as.integer(t),
+ as.integer(n),
+ as.integer(k)
+ )
+ m<-matrix(out$mat,nr=n)
+ # replace NAs with bitpattern corresponding to -2^31,
+ # i.e., 0 0 0.... 0 1
+ if (nas) {
+ idx1<-trunc(idx/nwords)+idx%%nwords # index for rows
+ idx2<-32+32*(1-idx%%2)*(nwords-1) # index for column
+ m[idx1,idx2]<-1
+ }
+ # calls user function to calculate statistic(s)
+ RET<-userfunc(m,...)
+ RET
+}
+
diff --git a/R/summary.RSctr.R b/R/summary.RSctr.R
new file mode 100755
index 0000000..c786819
--- /dev/null
+++ b/R/summary.RSctr.R
@@ -0,0 +1,13 @@
+"summary.RSctr" <-
+function(object,...)
+{
+ cat("\nCurrent sampler control specifications in",
+ deparse(substitute(object)),":\n")
+ cat("\tburn_in = ",object$burn_in,"\n")
+ cat("\tn_eff = ",object$n_eff,"\n")
+ cat("\tstep = ",object$step,"\n")
+ cat("\tseed = ",object$seed,"\n")
+ cat("\ttfixed = ",object$tfixed,"\n\n")
+ invisible(object)
+}
+
diff --git a/R/summary.RSmpl.R b/R/summary.RSmpl.R
new file mode 100755
index 0000000..e581c31
--- /dev/null
+++ b/R/summary.RSmpl.R
@@ -0,0 +1,16 @@
+"summary.RSmpl" <-
+function(object,...)
+{
+ cat("\nStatus of object",deparse(substitute(object)),
+ "after call to RSampler:\n")
+ cat("\tn = ",object$n,"\n")
+ cat("\tk = ",object$k,"\n")
+ cat("\tburn_in = ",object$burn_in,"\n")
+ cat("\tn_eff = ",object$n_eff,"\n")
+ cat("\tstep = ",object$step,"\n")
+ cat("\tseed = ",object$seed,"\n")
+ cat("\ttfixed = ",object$tfixed,"\n")
+ cat("\tn_tot = ",object$n_tot,"\n")
+ cat("\toutvec contains",length(object$outvec),"elements\n\n")
+}
+
diff --git a/R/summary.RSmplext.R b/R/summary.RSmplext.R
new file mode 100755
index 0000000..35f72e1
--- /dev/null
+++ b/R/summary.RSmplext.R
@@ -0,0 +1,15 @@
+"summary.RSmplext" <-
+function(object,...)
+{
+ cat("\nStatus of extracted object",
+ deparse(substitute(object)),":\n")
+ cat("\tn = ",object$n,"\n")
+ cat("\tk = ",object$k,"\n")
+ cat("\tburn_in = ",object$burn_in,"\n")
+ cat("\tn_eff = ",object$n_eff,"\n")
+ cat("\tstep = ",object$step,"\n")
+ cat("\tseed = ",object$seed,"\n")
+ cat("\ttfixed = ",object$tfixed,"\n")
+ cat("\tn_tot = ",object$n_tot,"\n")
+ cat("\toutvec contains",length(object$outvec),"elements\n\n")
+}
diff --git a/data/xmpl.rda b/data/xmpl.rda
new file mode 100755
index 0000000..ec3b6e9
Binary files /dev/null and b/data/xmpl.rda differ
diff --git a/data/xmplbig.rda b/data/xmplbig.rda
new file mode 100755
index 0000000..92d5d1b
Binary files /dev/null and b/data/xmplbig.rda differ
diff --git a/man/RSctr.Rd b/man/RSctr.Rd
new file mode 100755
index 0000000..d525254
--- /dev/null
+++ b/man/RSctr.Rd
@@ -0,0 +1,45 @@
+\name{RSctr}
+\alias{RSctr}
+\title{Control Object}
+\description{
+ The object of class \code{RSctr} represents the control parameter
+ specification for the sampling function \code{\link{rsampler}}.
+}
+\value{
+ A legitimate \code{\link{RSctr}} object is a list with components
+ \item{burn_in}{
+ the number of matrices to be sampled to
+ come close to a stationary distribution.
+ }
+ \item{n_eff}{
+ the number of effective matrices, i.e.,
+ the number of matrices
+ to be generated by the sampling function \code{\link{rsampler}}.
+ }
+ \item{step}{
+ controls the number number of void matrices generated in the the burn in
+ process and when effective matrices are generated (see note
+ in \code{\link{rsctrl}}). }
+ \item{seed}{
+ is the indicator for the seed of the random number generator.
+ If the value of seed at equals zero, a seed is generated
+ by the sampling function \code{\link{rsampler}}
+ }
+ \item{tfixed}{
+ \code{TRUE} or \code{FALSE}. \code{tfixed = TRUE} has no effect
+ if the input matrix is not quadratic,
+ i.e., all matrix elements are considered free (unrestricted).
+ If the input matrix is quadratic, and \code{tfixed = TRUE},
+ the main diagonal of the matrix is considered as fixed.
+ }
+}
+\section{Generation}{
+ This object is returned from function
+ \code{rsctrl}.
+}
+\section{Methods}{
+ This class has a method for the generic \code{summary}
+ function.
+}
+\seealso{\code{\link{rsctrl}} }
+\keyword{misc}
diff --git a/man/RSmpl.Rd b/man/RSmpl.Rd
new file mode 100755
index 0000000..2180709
--- /dev/null
+++ b/man/RSmpl.Rd
@@ -0,0 +1,46 @@
+\name{RSmpl}
+\alias{RSmpl}
+\alias{RSmplext}
+\title{Sample Objects}
+\description{
+ The objects of class \code{RSmpl} and \code{RSmplext} contain
+ the original input matrix, the generated (encoded) random matrices, and
+ some information about the sampling process.
+}
+\value{
+ A list of class \code{RSmpl} or \code{RSmplext} with components
+ \item{n}{number of rows of the input matrix}
+ \item{k}{number of columns of the input matrix}
+ \item{inpmat}{the input matrix}
+ \item{tfixed}{\code{TRUE}, if diagonals of \code{inpmat} are fixed}
+ \item{burn_in}{length of the burn in process}
+ \item{n_eff}{number of generated matrices (effective matrices)}
+ \item{step}{controls the number number of void matrices generated in the the burn in
+ process and when effective matrices are generated (see note
+ in \code{\link{rsctrl}}). }
+ \item{seed}{starting value for the random number generator}
+ \item{n_tot}{number of matrices in \code{outvec}.}
+ \item{outvec}{vector of encoded random matrices}
+ \item{ier}{error code (see below)}
+}
+\note{By default, all generated matrices plus
+ the original matrix (in position 1) are contained in
+ \code{outvec}, thus \code{n_tot = n_eff + 1}. If
+ the original matrix is not in \code{outvec} then
+ \code{n_tot = n_eff}.\cr\cr
+ If \code{ier} is 0, no error was detected. Otherwise use
+ the error function \code{rserror(ier)} to obtain some informations.\cr\cr
+ For saving and loading objects
+ of class \code{RSmpl} or \code{RSmplext}
+ see the example in \code{\link{rsextrobj}}.
+}
+\section{Generation}{
+ These classes of objects are returned from
+ \code{rsampler} and \code{rsextrobj}.
+}
+\section{Methods}{
+ Both classes have methods for the generic \code{summary}
+ function.
+}
+\seealso{\code{\link{rsampler}}, \code{\link{rsextrobj}} }
+\keyword{misc}
diff --git a/man/RaschSampler.package.Rd b/man/RaschSampler.package.Rd
new file mode 100755
index 0000000..1781455
--- /dev/null
+++ b/man/RaschSampler.package.Rd
@@ -0,0 +1,80 @@
+\name{RaschSampler-package}
+\alias{RaschSampler-package}
+\alias{RaschSampler}
+\docType{package}
+\title{Rasch Sampler Package}
+\description{
+The package implements an MCMC algorithm for sampling of
+binary matrices with fixed margins complying to the Rasch model.
+Its stationary distribution is uniform. The algorithm also allows
+for square matrices with fixed diagonal.\cr
+
+Parameter estimates in the Rasch model only depend on the marginal totals of
+the data matrix that is used for the estimation. From this it follows that, if the
+model is valid, all binary matrices with the same marginals as the observed one
+are equally likely. For any statistic of the data matrix, one can approximate
+the null distribution, i.e., the distribution if the Rasch model is valid, by taking
+a random sample from the collection of equally likely data matrices and constructing
+the observed distribution of the statistic.
+One can then simply determine the exceedence probability of the statistic in the
+observed sample, and thus construct a non-parametric test of the Rasch model.
+The main purpose of this package is the implementation of a methodology to build nonparametric
+tests for the Rasch model. \cr
+
+In the context of social network theories, where the structure of binary asymmetric
+relations is studied, for example,
+person \eqn{a} esteems person \eqn{b}, which correponds to a 1 in cell \eqn{(a, b)}
+of the associated adjacency matrix. If one wants to study
+the distribution of a statistic defined on the adjacency matrix and conditional
+on the marginal totals, one has to exclude the diagonal cells from consideration, i.e.,
+by keeping the diagonal cells fixed at an arbitrary value. The \code{RaschSampler} package
+has implemented an appropriate option, thus it can be also used for sampling random adjacency
+matrices with given marginal totals.
+}
+\details{
+\tabular{ll}{
+Package: \tab RaschSampler\cr
+Type: \tab Package\cr
+Version: \tab 0.8-4\cr
+Date: \tab 2010-03-22\cr
+License: \tab GNU GPL 2, June 1991\cr
+}
+The user has to supply a binary input matrix. After defining appropriate control
+parameters using \code{\link{rsctrl}} the sampling function \code{\link{rsampler}}
+may be called to obtain an object of class \code{\link{RSmpl}} which contains the
+generated random matrices in encoded form. After defining an appropriate function
+to operate on a binary matrix (e.g., calculate a statistic such as \code{\link{phi.range}})
+the application of this function to the sampled matrices is performed
+using \code{\link{rstats}}. Prior to applying the user defined function, \code{\link{rstats}}
+decodes the matrices packed in the \code{\link{RSmpl}}-object.\cr
+
+The package also defines a utility function \code{\link{rsextrobj}} for extracting certains parts from
+the \code{\link{RSmpl}}-object resulting in an object of class \code{\link{RSmplext}}.
+Both types of objects can be saved and reloaded for later use.\cr
+
+Summary methods are available to print information on these objects, as well as
+on the control object \code{\link{RSctr}} which is obtained from using
+\code{\link{rsctrl}} containing the specification for the sampling routine.
+
+}
+\author{Reinhold Hatzinger, Patrick Mair, Norman Verhelst
+
+Maintainer: <reinhold.hatzinger at wu-wien.ac.at>
+}
+\references{
+Verhelst, N. D. (2008) An Efficient MCMC Algorithm to Sample Binary
+Matrices with Fixed Marginals. Psychometrika, Volume 73, Number 4\cr
+Verhelst, N. D., Hatzinger, R., and Mair, P. (2007) The Rasch Sampler, Journal of Statistical Software, Vol. 20, Issue 4, Feb 2007
+}
+\note{The current implementation allows for data matrices up to 4096 rows and 128 columns.
+ This can be changed by setting \code{nmax} and \code{kmax} in \code{RaschSampler.f90}
+ to values which are a power of 2. These values should also be changed in \code{rserror.R}.
+
+ For convenience, we reuse the Fortran code of package version 0.8-1 which cicumvents the
+ compiler bug in Linux distributions of GCC 4.3. The following note from package version 0.8-3
+ is thus obsolete:
+ In case of compilation errors (due to a bug in Linux distributions of GCC 4.3) please use
+ \code{RaschSampler.f90} from package version 0.8-1 and change \code{nmax} and \code{kmax}
+ accordingly (or use GCC 4.4).}
+\keyword{ package }
+
diff --git a/man/phi.range.Rd b/man/phi.range.Rd
new file mode 100755
index 0000000..38436e4
--- /dev/null
+++ b/man/phi.range.Rd
@@ -0,0 +1,25 @@
+\name{phi.range}
+\alias{phi.range}
+\title{ Example User Function }
+\description{
+ Calculates the \eqn{R_\phi} statistic, i.e., the range of
+ the inter-column correlations (\eqn{\phi}-coefficients) for
+ a binary matrix.
+}
+\usage{
+phi.range(mat)
+}
+\arguments{
+ \item{mat}{ a binary matrix }
+}
+\value{
+ the range of the inter-column correlations
+}
+\examples{
+ctr <- rsctrl(burn_in = 10, n_eff = 5, step=10, seed = 123, tfixed = FALSE)
+mat <- matrix(sample(c(0,1), 50, replace = TRUE), nr = 10)
+rso <- rsampler(mat, ctr)
+rso_st <- rstats(rso,phi.range)
+print(unlist(rso_st))
+}
+\keyword{misc}
diff --git a/man/rsampler.Rd b/man/rsampler.Rd
new file mode 100755
index 0000000..c59d48d
--- /dev/null
+++ b/man/rsampler.Rd
@@ -0,0 +1,63 @@
+\name{rsampler}
+\alias{rsampler}
+\title{Sampling Binary Matrices}
+\description{
+ The function implements an MCMC algorithm for sampling of binary
+ matrices with fixed margins complying to the Rasch model.
+ Its stationary distribution is
+ uniform. The algorithm also allows for square matrices with
+ fixed diagonal.
+}
+\usage{
+rsampler(inpmat, controls = rsctrl())
+}
+\arguments{
+ \item{inpmat}{ A binary (data) matrix with \eqn{n} rows and \eqn{k} columns.}
+ \item{controls}{An object of class \code{\link{RSctr}}. If not specified, the default
+ parameters as returned by function \code{\link{rsctrl}} are
+ used.}
+}
+\details{
+ \code{rsampler} is a wrapper function for a Fortran routine to generate binary random matrices based
+ on an input matrix.
+ On output the generated binary matrices are integer encoded. For further
+ processing of the generated matrices use the function \code{\link{rstats}}.
+}
+\value{
+ A list of class \code{\link{RSmpl}} with components
+ \item{n}{number of rows of the input matrix}
+ \item{k}{number of columns of the input matrix}
+ \item{inpmat}{the input matrix}
+ \item{tfixed}{\code{TRUE}, if diagonals of \code{inpmat} are fixed}
+ \item{burn_in}{length of the burn in process}
+ \item{n_eff}{number of generated matrices (effective matrices)}
+ \item{step}{controls the number number of void matrices generated in the the burn in
+ process and when effective matrices are generated (see note
+ in \code{\link{rsctrl}}). }
+ \item{seed}{starting value for the random number generator}
+ \item{n_tot}{number of matrices in \code{outvec}, \code{n_tot = n_eff + 1}}
+ \item{outvec}{vector of encoded random matrices}
+ \item{ier}{error code}
+}
+\references{Verhelst, N. D. (2008) An Efficient MCMC Algorithm to Sample Binary
+Matrices with Fixed Marginals. Psychometrika, Volume 73, Number 4}
+\author{Reinhold Hatzinger, Norman Verhelst}
+\note{
+ An element of \code{outvec} is a four byte (or 32 bits) integer. The matrices
+ to be output are stored bitwise (some bits
+ are unused, since a integer is used for every row of a matrix. So
+ the number of integers per row needed equals (k+31)/32 (integer division),
+ which is one to four in the present implementation since the number of columns
+ and rows must not exceed 128 and 4096, respectively.\cr
+
+ The summary method (\code{\link{summary.RSmpl}}) prints
+ information on the content of the output object.
+}
+\seealso{\code{\link{rsctrl}}, \code{\link{rstats}} }
+\examples{
+data(xmpl)
+ctr<-rsctrl(burn_in=10, n_eff=5, step=10, seed=0, tfixed=FALSE)
+res<-rsampler(xmpl,ctr)
+summary(res)
+}
+\keyword{misc}
diff --git a/man/rsctrl.Rd b/man/rsctrl.Rd
new file mode 100755
index 0000000..ccf1b11
--- /dev/null
+++ b/man/rsctrl.Rd
@@ -0,0 +1,86 @@
+\name{rsctrl}
+\alias{rsctrl}
+\title{Controls for the Sampling Function}
+\description{
+
+ Various parameters that control aspects of
+ the random generation of binary matrices.
+}
+\usage{
+rsctrl(burn_in = 100, n_eff = 100, step = 16, seed = 0, tfixed = FALSE)
+}
+\arguments{
+ \item{burn_in}{
+ the number of sampled matrices to
+ come close to a stationary distribution.
+ The default is \code{burn_in = 100}.
+ (The actual number is \code{2 * burn_in * step}.)
+ }
+ \item{n_eff}{
+ the number of effective matrices, i.e.,
+ the number of matrices
+ to be generated by the sampling function \code{\link{rsampler}}.
+ \code{n_eff} must be positive and not larger than 8191
+ (2\eqn{\mbox{\textasciicircum}}{^}13-1).
+ The default is \code{n_eff = 100}.
+ }
+ \item{step}{controls the number number of void matrices generated in the the burn in
+ process and when effective matrices are generated (see note
+ below). The default is \code{step = 16}. }
+ \item{seed}{
+ is the indicator for the seed of the random number generator.
+ Its value must be in the range 0 and 2147483646 (2**31-2).
+ If the value of seed equals zero, a seed is generated
+ by the sampling function \code{\link{rsampler}}
+ (dependent on the system's clock) and its value is returned
+ in the output. If seed is not equal to zero, its
+ value is used as the seed of the random number generator.
+ In that case its value is unaltered at output.
+ The default is \code{seed = 0}.
+ }
+ \item{tfixed}{logical, -- specifies if in case of a quadratic input
+ matrix the diagonal is considered fixed (see note below).
+ The default is \code{tfixed = FALSE}.
+ }
+}
+
+\value{
+ A list of class \code{RSctr} with components
+ \code{burn_in}, \code{n_eff}, \code{step},
+ \code{seed}, \code{tfixed}.,
+}
+\note{
+ If one of the components is incorrectly specified
+ the error function \code{rserror}
+ is called and some informations are printed. The ouput object
+ will not be defined.\cr\cr
+ The specification of \code{step} controls the sampling algorithm as follows:
+ If , e.g., \code{burn_in = 10}, \code{n_eff = 5}, and \code{step = 2},
+ then during the burn in period \code{step * burn_in = 2 * 10}
+ matrices are generated. After that, \code{n_eff * step = 5 * 2} matrices
+ are generated and every second matrix of these last ten is returned from
+ \code{link{rsampler}}.\cr\cr
+ \code{tfixed} has no effect if the input matrix is not quadratic,
+ i.e., all matrix elements are considered free (unrestricted).
+ If the input matrix is quadratic, and \code{tfixed = TRUE},
+ the main diagonal of the matrix is considered as fixed.
+ On return from \code{link{rsampler}} all diagonal elements
+ of the generated matrices are set to zero.
+ This specification applies, e.g.,
+ to analyzing square incidence matrices
+ representing binary asymmetric relation
+ in social network theory.\cr\cr
+ The summary method (\code{\link{summary.RSctr}}) prints
+ the current definitions. \cr
+}
+\seealso{\code{\link{rsampler}} }
+\examples{
+ctr <- rsctrl(n_eff = 1, seed = 987654321) # specify new controls
+summary(ctr)
+
+\dontrun{
+ctr2 <- rsctrl(step = -3, n_eff = 10000) # incorrect specifications
+}
+}
+
+\keyword{misc}
diff --git a/man/rsextrmat.Rd b/man/rsextrmat.Rd
new file mode 100755
index 0000000..e97e553
--- /dev/null
+++ b/man/rsextrmat.Rd
@@ -0,0 +1,28 @@
+\name{rsextrmat}
+\alias{rsextrmat}
+\title{Extracting a Matrix}
+\description{
+ Convenience function to extract a matrix.
+}
+\usage{
+rsextrmat(RSobj, mat.no = 1)
+}
+\arguments{
+ \item{RSobj}{object as obtained from using \code{rsampler} or \code{rsextrobj}}
+ \item{mat.no}{number of the matrix to extract from the sample object.}
+}
+\value{
+ One of the matrices (either the original or a sampled matrix)
+}
+\seealso{\code{\link{rsampler}}, \code{\link{rsextrobj}},\code{\link{rstats}},}
+\examples{
+ctr <- rsctrl(burn_in = 10, n_eff = 3, step=10, seed = 0, tfixed = FALSE)
+mat <- matrix(sample(c(0,1), 50, replace = TRUE), nr = 10)
+all_m <- rsampler(mat, ctr)
+summary(all_m)
+
+# extract the third sampled matrix (here the fourth)
+third_m <- rsextrmat(all_m, 4)
+head(third_m)
+}
+\keyword{misc}
diff --git a/man/rsextrobj.Rd b/man/rsextrobj.Rd
new file mode 100755
index 0000000..8cef36f
--- /dev/null
+++ b/man/rsextrobj.Rd
@@ -0,0 +1,67 @@
+\name{rsextrobj}
+\alias{rsextrobj}
+\title{Extracting Encoded Sample Matrices}
+\description{
+ Utility function to extract some of the generated matrices, still in encoded form.
+}
+\usage{
+rsextrobj(RSobj, start = 1, end = 8192)
+}
+\arguments{
+ \item{RSobj}{object as obtained from using \code{rsampler}}
+ \item{start}{number of the matrix to start with. When specifying 1
+ (the default value) the original input matrix is
+ included in the output object.
+ }
+ \item{end}{last matrix to be extracted. If \code{end}
+ is not specified, all matrices from \code{RSobj}
+ are extracted (the maximal value is 8192, see
+ \code{rsctrl}). If \code{end} is larger than
+ the number of matrices stored in \code{RSobj},
+ \code{end} is set to the highest possible value
+ (i.e., \code{n_tot}).
+ }
+}
+\value{
+ A list of class \code{\link{RSmpl}} with components
+ \item{n}{number of rows of the input matrix}
+ \item{k}{number of columns of the input matrix}
+ \item{inpmat}{the input matrix}
+ \item{tfixed}{\code{TRUE}, if diagonals of \code{inpmat} are fixed}
+ \item{burn_in}{length of the burn in process}
+ \item{n_eff}{number of generated matrices (effective matrices)}
+ \item{step}{controls the number number of void matrices generated in the burn in
+ process and when effective matrices are generated (see note
+ in \code{\link{rsctrl}}). }
+ \item{seed}{starting value for the random number generator}
+ \item{n_tot}{number of matrices in \code{outvec}.}
+ \item{outvec}{vector of encoded random matrices}
+ \item{ier}{error code}
+}
+\note{By default, all generated matrices plus
+ the original matrix (in position 1) are contained in
+ \code{outvec}, thus \code{n_tot = n_eff + 1}. If
+ the original matrix is not in \code{outvec} then
+ \code{n_tot = n_eff}.\cr
+ For saving and loading objects
+ of class \code{RSobj} see the example below.
+
+ For extracting a decoded (directly usable) matrix use \code{\link{rsextrmat}}.
+}
+\seealso{\code{\link{rsampler}}, \code{\link{rsextrmat}} }
+\examples{
+ctr <- rsctrl(burn_in = 10, n_eff = 3, step=10, seed = 0, tfixed = FALSE)
+mat <- matrix(sample(c(0,1), 50, replace = TRUE), nr = 10)
+all_m <- rsampler(mat, ctr)
+summary(all_m)
+
+some_m <- rsextrobj(all_m, 1, 2)
+summary(some_m)
+
+\dontrun{
+save(some_m, file = "some.RSobj")
+some_new <- load("some.RSobj")
+summary(some_new)
+}
+}
+\keyword{misc}
diff --git a/man/rstats.Rd b/man/rstats.Rd
new file mode 100755
index 0000000..5bc6807
--- /dev/null
+++ b/man/rstats.Rd
@@ -0,0 +1,66 @@
+\name{rstats}
+\alias{rstats}
+\title{Calculating Statistics for the Sampled Matrices}
+\description{
+ This function is used to calculate user defined statistics for the
+ (original and) sampled matrices. A user defined function has to
+ be provided.
+}
+\usage{
+rstats(RSobj, userfunc, ...)
+}
+\arguments{
+ \item{RSobj}{object as obtained from using \code{\link{rsampler}}
+ or \code{\link{rsextrobj}} }
+ \item{userfunc}{a user defined function which performs operations
+ on the (original and) sampled matrices. The first argument in the definition
+ of the user function must be an object of type matrix.}
+ \item{...}{further arguments, that are passed to the user function}
+}
+\value{
+ A list of objects as specified in the user supplied function
+}
+\note{The encoded matrices that are contained in the
+ input object \code{RSobj} are decoded and passed to the user function in turn.
+ If \code{RSobj} is not an object obtained from either \code{\link{rsampler}}
+ or \code{\link{rsextrobj}} or
+ no user function is specified an error message is printed.
+ A simple user function, \code{\link{phi.range}}, is included in
+ the RaschSampler package for demonstration purposes.\cr
+
+ \code{rstats} can be used to obtain the 0/1 values for any
+ of the sampled matrices (see second example below). Please note,
+ that the output from the user function is stored in a list where
+ the number of components corresponds to the number of matrices passed
+ to the user function (see third example).
+}
+\seealso{\code{\link{rsampler}}, \code{\link{rsextrobj}} }
+\examples{
+ctr <- rsctrl(burn_in = 10, n_eff = 5, step=10, seed = 12345678, tfixed = FALSE)
+mat <- matrix(sample(c(0,1), 50, replace = TRUE), nr = 10)
+rso <- rsampler(mat, ctr)
+rso_st <- rstats(rso,phi.range)
+unlist(rso_st)
+
+# extract the third generated matrix
+# (here, the first is the input matrix)
+# and decode it into rsmat
+
+rso2 <- rsextrobj(rso,4,4)
+summary(rso2)
+rsmat <- rstats(rso2, function(x) matrix(x, nr = rso2$n))
+print(rsmat[[1]])
+
+# extract only the first r rows of the third generated matrix
+
+mat<-function(x, nr = nr, r = 3){
+ m <- matrix(x, nr = nr)
+ m[1:r,]
+}
+rsmat2 <- rstats(rso2, mat, nr=rso$n, r = 3)
+print(rsmat2[[1]])
+
+# apply a user function to the decoded object
+print(phi.range(rsmat[[1]]))
+}
+\keyword{misc}
diff --git a/man/summary.RSctr.Rd b/man/summary.RSctr.Rd
new file mode 100755
index 0000000..5d0a0a8
--- /dev/null
+++ b/man/summary.RSctr.Rd
@@ -0,0 +1,19 @@
+\name{summary.RSctr}
+\alias{summary.RSctr}
+\title{Summary Method for Control Objects}
+\description{
+ Prints the current definitions for the sampling function.
+}
+\usage{
+\method{summary}{RSctr}(object, ...)
+}
+\arguments{
+ \item{object}{ object of class \code{RSctr} as obtained from \code{\link{rsctrl}} }
+ \item{\dots}{ potential further arguments (ignored) }
+}
+\seealso{ \code{\link{rsctrl}} }
+\examples{
+ ctr <- rsctrl(n_eff = 1, seed = 123123123) # specify controls
+ summary(ctr)
+}
+\keyword{misc}
diff --git a/man/summary.RSmpl.Rd b/man/summary.RSmpl.Rd
new file mode 100755
index 0000000..5df2509
--- /dev/null
+++ b/man/summary.RSmpl.Rd
@@ -0,0 +1,31 @@
+\name{summary.RSmpl}
+\alias{summary.RSmpl}
+\alias{summary.RSmplext}
+\title{Summary Methods for Sample Objects}
+
+\description{
+ Prints a summary list for sample objects of class \code{\link{RSmpl}}
+ and \code{\link{RSmplext}}.
+}
+\usage{
+\method{summary}{RSmpl}(object, ...)
+\method{summary}{RSmplext}(object, ...)
+}
+\arguments{
+ \item{object}{object as obtained from \code{rsampler} or \code{rsextrobj} }
+ \item{\dots}{ potential further arguments (ignored) }
+}
+\details{
+ Describes the status of an sample object.
+}
+\seealso{\code{\link{rsampler}}, \code{\link{rsextrobj}} }
+\examples{
+ctr <- rsctrl(burn_in = 10, n_eff = 3, step=10, seed = 0, tfixed = FALSE)
+mat <- matrix(sample(c(0,1), 50, replace = TRUE), nr = 10)
+all_m <- rsampler(mat, ctr)
+summary(all_m)
+
+some_m <- rsextrobj(all_m, 1, 2)
+summary(some_m)
+}
+\keyword{misc}
diff --git a/man/xmpl.Rd b/man/xmpl.Rd
new file mode 100755
index 0000000..218309b
--- /dev/null
+++ b/man/xmpl.Rd
@@ -0,0 +1,25 @@
+\name{xmpl}
+\alias{xmpl}
+\alias{xmplbig}
+\docType{data}
+\title{Example Data}
+\description{
+ Ficitious data sets - matrices with binary responses
+}
+\usage{data(xmpl)}
+\format{
+ The format of \code{xmpl} is:\cr
+ 300 rows (referring to subjects) \cr
+ 30 columns (referring to items) \cr
+
+ The format of \code{xmplbig} is:\cr
+ 4096 rows (referring to subjects) \cr
+ 128 columns (referring to items) \cr
+ \code{xmplbig} has the maximum dimensions that the RaschSampler package
+ can handle currently.
+}
+\examples{
+data(xmpl)
+print(head(xmpl))
+}
+\keyword{datasets}
diff --git a/src/RaschSampler.f90 b/src/RaschSampler.f90
new file mode 100755
index 0000000..0547815
--- /dev/null
+++ b/src/RaschSampler.f90
@@ -0,0 +1,550 @@
+subroutine sampler(n,k,inputmat,tfixed,burn_in,n_eff,step,seed,outputvec,ier)
+
+ ! sample binary matrices with given marginals.
+ ! input: n = number of rows (integer*4)
+ ! k = number of columns (integer*4)
+ ! inputmat: input binary matrix (n*k) (integer*4)
+ ! tfixed: main diagonal is fixed if true; has only effect if n.eq.k
+ ! step: if the matrix #i is an effective matrix, then the next effective matrix is #(i+step)
+ ! burn_in: number of burn_in matrices in units of step
+ ! n_eff: number of effective matrices, to be written in the output vector
+ ! I/O seed (integer*4). Seed of the random generator
+ ! if seed <> 0: seed is unaltered and |seed| is used as the seed of the random generator.
+ ! if seed.eq.0: a seed is generated from the system clock, and its value is returned.
+ ! ATTENTION: currently seed.eq.0 is deactivated and only seed<>0 is allowed
+ ! (see lines 103-113) (rh 2006-10-25)
+ ! output:outputvec (integer*4 vector): n_eff binary matrices, stored bitwise in the following way:
+ ! if(k<=32) one row of a matrix is stored in one position (number) of four bytes, the first element
+ ! is bit 0, the second is bit 1, etc.
+ ! if(32<k<=64) two positions are used per row: the frst 32 in position 1, the remainder in position 2,
+ ! again starting with element 33 in bit 0, element 34 in bit 1, etc.
+ ! not used bits are set to zero.
+ ! ier: output error code
+ ! ier = 0: O.K.
+ ! 1: n > nmax = 1024 = 2**10 !!!!!! changed to 2**12 in 0.8-3
+ ! 2: k > kmax = 64 = 2**6 !!!!!! changed to 2**7 in 0.8-3
+ ! 4: n_eff > n_effmax = 8191 = 2**13 - 1
+ ! 8: burn_in < 0
+ ! 16: step <= 0
+ ! 1-31: sums of the foregoing codes
+ ! 32: input matrix contains values other than one or zero
+ ! 64: the input matrix has a Guttman form
+ ! if tfixed and n.eq.k, the main diagonal is condidered as fixed.
+ ! if tfixed and n!=k, the case is treated as a rectangular matrix without constraints
+ ! the Markov chain transition matrix used is Q**step
+
+ integer(kind=4), dimension(n*(k+31)/32*(n_eff+1)),intent(out)::outputvec
+ integer(kind=4), intent(out) :: ier
+ integer(kind=4), intent(in) :: n,k,burn_in,n_eff,step
+ integer(kind=4), dimension(n,k),intent(in) :: inputmat
+ integer(kind=4) :: seed
+ logical(kind=4), intent(in) :: tfixed
+
+ character(len=10) :: timevec
+
+ !integer(kind=4),parameter :: nmax=1024,kmax=64,n_effmax=8191 !!!!!! kmax changed to 2**7 nmax changed to 2**12
+ integer(kind=4),parameter :: nmax=4096,kmax=128,n_effmax=8191
+ integer(kind=4), allocatable :: a(:),b(:),aold(:),bold(:),a_kol(:),b_kol(:),iwork(:)
+ integer(kind=4) :: i,j,m,kk2,kk3,it,krand,k2,k3,nhex,k2old,k3old,nhexold
+ integer(kind=4) :: x1,x2 ! x1 and x2 are reserved for the random generator
+ integer(kind=4) :: words_per_row, offset
+ integer(kind=4),allocatable :: hexmat(:,:),hexmatold(:,:)
+
+ real(kind=4) :: tijd
+
+ logical(kind=1),parameter :: t=.true.,f=.false.
+ logical(kind=1),allocatable :: t_in(:,:)
+ logical(kind=1),dimension(3,3) :: hexa,hexb
+ logical(kind=1),allocatable :: twa(:),twb(:),tw(:),tng(:),tngold(:),col1(:),col2(:) !tng = non-guttman pair
+ logical(kind=1) :: t_eff,tfixnow
+
+ data hexa/f,f,t,t,f,f,f,t,f/,hexb/f,t,f,f,f,t,t,f,f/
+
+ !check error codes 1, 2 4, 8 and 16
+ ier=0
+ if(n.le.0 .or. n.gt.nmax)ier=ier+1
+ if(k.le.0 .or. k.gt.kmax)ier=ier+2
+ if(n_eff.le.0 .or. n_eff.gt.n_effmax)ier=ier+4
+ if(burn_in.lt.0)ier=ier+8
+ if(step.le.0)ier=ier+16
+ if(ier.ne.0)return
+
+ ! allocate the necessary arrays
+
+ kk2=k*(k-1)/2
+ allocate(t_in(n,k))
+ allocate(a(kk2),b(kk2),aold(kk2),bold(kk2),a_kol(kk2),b_kol(kk2))
+ allocate(twa(n),twb(n),tw(n),tng(kk2),tngold(kk2),col1(n),col2(n))
+ allocate (iwork(n))
+
+ tfixnow=tfixed
+ if(n.ne.k)tfixnow=.false.
+ if(tfixnow)then
+ kk3=kk2*(k-2)/3
+ allocate(hexmat(3,kk3),hexmatold(3,kk3))
+ endif
+
+ ! check error code 32
+ ier=count(inputmat.gt.1)+count(inputmat.lt.0)
+ if(ier.ne.0)then
+ ier=32
+ return
+ endif
+
+ ! copy input matrix to t_in
+ !!!t_in=inputmat
+ !!!replaced by
+ t_in=btest(inputmat,0)
+ if(tfixnow) then
+ forall (i=1:n) t_in(i,i)=f
+ endif
+! !select seed for random number generation
+! if(seed.eq.0)then
+! call date_and_time(TIME=timevec)
+! read(timevec,'(f10.3)')tijd
+! x1=tijd*1000.
+! x1=0.
+! x1=x1+536870911 ! =x1 + 2**29 - 1
+! seed=x1
+! else
+! x1=abs(seed)
+! endif
+ x1=abs(seed) ! added from upper else clause (to be removed if random seed enabled)
+ x2=x1
+ call rand1(x1)
+ krand=1+mod(k,25) ! KRAND selects a random generator (used in RAND_INTEGER42)
+
+ !fill the arrays a_kol, b_kol, a, b, and b_pairs for the input matrices t_in
+ !determine the weight k2 (= #neighbour column pairs of the input matrix)
+ it=0
+ do i=2,k
+ do j=1,i-1
+ it=it+1
+ a_kol(it)=i
+ b_kol(it)=j
+ call findab(t_in(1:n,i),t_in(1:n,j),i,j,a(it),b(it))
+ tng(it)=a(it)*b(it).gt.0
+ end do
+ end do
+ k2=count(tng(1:kk2))
+ k3=k2
+ if(tfixnow)then
+ call hexagon
+ if(nhex.gt.0)k3=k2+1
+ endif
+ ! check on the Guttman condition (error code 64)
+ if(k3.eq.0)then
+ ier=64
+ return
+ endif
+
+ ! pack the input matrix and put it in the outputvector
+ offset=0 ! offset is the number of words defined until now
+ words_per_row=(k+31)/32
+ call pack_matrix(outputvec) ! offset is updated within the subroutine
+
+ ! compute the total number of matrices to be computed
+ n_tot=burn_in+n_eff
+ ! do the sampling
+ do i=1,n_tot
+ t_eff=i.gt.burn_in
+ do j=1,step ! generate step matrices before computing the statistic
+ call rand_integer42(it,k3,x1,x2,krand)
+ !{*** from here to ***} only applies to alternating hexagons
+ if(k3.gt.k2.and.it.eq.k3)then ! there are restrictions on the main diagonal
+ ! and there exists at least one alternating hexagon (k3>k2)
+ ! an alternating hexagon has to be changed (it=k3)
+ ! first save necessary elements for Metropolis-Hastings
+ k2old=k2
+ k3old=k3
+ aold=a
+ bold=b
+ nhexold=nhex
+ hexmatold(:,1:nhex)=hexmat(:,1:nhex)
+ tngold=tng
+ ! make a new matrix by switching one alternating hexagon
+ call rand_integer42(it,nhex,x1,x2,krand)
+ call make_matrix3(it)
+ call update_pairs3(hexmat(1,it),hexmat(2,it),hexmat(3,it))
+ call update_pairs3(hexmat(2,it),hexmat(1,it),hexmat(3,it))
+ call update_pairs3(hexmat(3,it),hexmat(1,it),hexmat(2,it))
+ call hexagon
+ k2=count(tng(1:kk2))
+ k3=k2+1 ! there is at least one alternating hexagon, vz, the one that just changed into its complement
+ if(k3old.lt.k3) then ! check if process possibly remains in the same state (Metropolis-Hastings)
+ call rand_integer42(m,k3,x1,x2,krand)
+ if(m.gt.k3old)then ! process remains in the same state
+ call make_matrix3(it) ! this restores the matrix
+ k3=k3old
+ k2=k2old
+ a=aold
+ b=bold
+ nhex=nhexold
+ hexmat(:,1:nhex)=hexmatold(:,1:nhex)
+ tng=tngold
+ endif
+ endif
+ cycle
+ endif
+ !***}
+ call pick_a_pair(it)
+ col1=t_in(1:n,a_kol(it))
+ col2=t_in(1:n,b_kol(it))
+ k2old=k2
+ k3old=k3
+ aold=a
+ bold=b
+ tngold=tng
+ if(tfixnow)then
+ nhexold=nhex
+ hexmatold(:,1:nhex)=hexmat(:,1:nhex)
+ endif
+ call make_matrix(it)
+ call update_pairs(a_kol(it),b_kol(it))
+ call update_pairs(b_kol(it),a_kol(it))
+ k2=count(tng(1:kk2))
+ k3=k2
+ if(tfixnow)then
+ call hexagon
+ if(nhex.gt.0)k3=k3+1
+ endif
+ if(k3old.lt.k3)then ! apply Metropolis-Hastings
+ call rand_integer42(m,k3,x1,x2,krand)
+ if(m.gt.k3old)then ! process remains in the current state
+ t_in(1:n,a_kol(it))=col1(1:n)
+ t_in(1:n,b_kol(it))=col2(1:n)
+ k2=k2old
+ k3=k3old
+ a=aold
+ b=bold
+ tng=tngold
+ if(tfixnow)then
+ nhex=nhexold
+ hexmat(:,1:nhex)=hexmatold(:,1:nhex)
+ endif
+ endif
+ endif
+ end do
+ ! if this is an 'effective' sample, pack it and store it in outputvec
+ if(t_eff)call pack_matrix(outputvec)
+ end do ! here ends the sampling procedure
+
+ deallocate(a,b,aold,bold,a_kol,b_kol,iwork,twa,twb,tw,tng,tngold,col1,col2,t_in)
+ if(tfixnow) deallocate(hexmat,hexmatold)
+
+
+ contains
+
+ subroutine pack_matrix(vec)
+ integer(kind=4) vec(*)
+ integer(kind=4) :: i,j,ib,ie,it,iw
+ do i=1,n
+ ib=1
+ do iw=1,words_per_row
+ offset=offset+1
+ vec(offset)=0
+ ie=min(ib+31,k)
+ it=-1
+ do j=ib,ie
+ it=it+1
+ if(.not.t_in(i,j))cycle
+ vec(offset)=ibset(vec(offset),it)
+ end do
+ ib=ie+1
+ end do
+ end do
+ end subroutine pack_matrix
+
+ subroutine findab(ta,tb,i,j,a,b)
+ !!!!! logical(kind=1):: ta(n),tb(n)
+ logical(kind=1):: ta(n),tb(n),test(n)
+ integer(kind=4)::a,b,i,j
+ tw=(ta.neqv.tb)
+ if(tfixnow)then
+ tw(i)=.false.
+ tw(j)=.false.
+ endif
+ test=ta.and.tw
+ a=count(test)
+ test=tb.and.tw
+ b=count(test)
+ !!!!! a=count(ta.and.tw)
+ !!!!! b=count(tb.and.tw)
+ end subroutine findab
+
+
+ subroutine pick_a_pair(it)
+ integer(kind=4),intent(out)::it
+ integer(kind=4) ::i,m
+ call rand_integer42(it,k2,x1,x2,krand)
+ m=count(tng(1:it))
+ if(m.eq.it)return
+ do i=it+1,kk2
+ if(.not.tng(i))cycle
+ m=m+1
+ if(m.eq.it)then
+ it=i
+ return
+ endif
+ end do
+ end subroutine pick_a_pair
+
+ subroutine make_matrix(it)
+ integer(kind=4),intent(in)::it
+ integer(kind=4) ::m,i,j,ii,jj
+!!!!! logical(kind=1),allocatable :: test(:)
+ logical(kind=1) :: test(n)
+
+ ii=a_kol(it)
+ jj=b_kol(it)
+ if(a(it)*b(it).eq.1)then ! columns ii and jj contain a single tetrad.
+ ! no sampling is necessary: the tetrad is complemented
+ j=0
+ do i=1,n
+ if(tfixnow)then
+ if(i.eq.ii)cycle
+ if(i.eq.jj)cycle
+ endif
+ if(t_in(i,ii).eqv.t_in(i,jj))cycle
+ j=j+1
+ t_in(i,ii)=t_in(i,jj)
+ t_in(i,jj)=.not.t_in(i,ii)
+ if(j.eq.2)return
+ end do
+ endif
+ do ! a random binomial operation is applied
+ ! copy the two selected colums in the logical vectors twa and twb
+ twa(1:n)=t_in(1:n,ii)
+ twb(1:n)=t_in(1:n,jj)
+ m=a(it)+b(it)
+ call combine(m,a(it),iwork,tw) ! generate the random combination of a(it) objets out of m
+ ! insert the combination into the vectors twa and twb
+ j=0
+ do i=1,n
+ if(tfixnow)then
+ if(i.eq.ii)cycle
+ if(i.eq.jj)cycle
+ endif
+ if(twa(i).eqv.twb(i))cycle
+ j=j+1
+ if(tw(j))then
+ twa(i)=.true.
+ twb(i)=.false.
+ else
+ twa(i)=.false.
+ twb(i)=.true.
+ endif
+ if(j.eq.m)exit
+ end do
+ ! check whether matrix has changed
+ !!!!!
+ !!!!! test=twa
+ test=twa(1:n).eqv.t_in(1:n,ii)
+ m=count(test)
+ !!!!! m=count(twa(1:n).eqv.t_in(1:n,ii))
+ if(m.ne.n)exit ! the matrix has changed
+ end do ! the matrix has not changed; a new combination is tried.
+ ! the changes are inserted in the matrix t_in
+ t_in(1:n,ii)=twa(1:n)
+ t_in(1:n,jj)=twb(1:n)
+
+ end subroutine make_matrix
+
+ subroutine make_matrix3(it)
+ integer(kind=4), intent(in)::it
+ integer(kind=4) ::i,j,ii,jj
+ do i=1,2
+ ii=hexmat(i,it)
+ do j=i+1,3
+ jj=hexmat(j,it)
+ t_in(ii,jj)=.not.t_in(ii,jj)
+ t_in(jj,ii)=.not.t_in(jj,ii)
+ end do
+ end do
+ end subroutine make_matrix3
+
+ subroutine combine(n,k,ix,tx)
+ ! generate a random combination of k objects out of n
+ ! the result is stored in the logical n-vector tx
+ ! ix is a working array
+ integer(kind=4) ::n,k,kk,ii,iu,nnu
+ integer(kind=4) ::ix(n)
+ logical(kind=1) ::tx(n)
+
+ ix(1:n)=(/(ii,ii=1,n)/)
+ tx(1:n)=.false.
+ nnu=n
+ kk=min(k,n-k)
+ do ii=1,kk
+ call rand_integer42(iu,nnu,x1,x2,krand)
+ tx(ix(iu))=.true.
+ if(iu.lt.nnu)ix(iu:nnu-1)=ix(iu+1:nnu)
+ nnu=nnu-1
+ end do
+ if(kk.lt.k)tx=.not.tx
+ end subroutine combine
+
+ subroutine update_pairs(i,j)
+ integer(kind=4),intent(in) :: i,j
+ integer(kind=4):: jt,m
+ do m=1,i-1
+ if(m.eq.j)cycle
+ jt=(i-1)*(i-2)/2+m
+ call findab(t_in(1:n,i),t_in(1:n,m),i,m,a(jt),b(jt))
+ tng(jt)=a(jt)*b(jt).gt.0
+ end do
+ do m=i+1,k
+ if(m.eq.j)cycle
+ jt=(m-1)*(m-2)/2+i
+ call findab(t_in(1:n,m),t_in(1:n,i),m,i,a(jt),b(jt))
+ tng(jt)=a(jt)*b(jt).gt.0
+ end do
+ end subroutine update_pairs
+
+ subroutine update_pairs3(i,j,l)
+ integer(kind=4),intent(in) :: i,j,l
+ integer(kind=4):: jt,m
+ do m=1,i-1
+ if(m.eq.j)cycle
+ if(m.eq.l)cycle
+ jt=(i-1)*(i-2)/2+m
+ call findab(t_in(1:n,i),t_in(1:n,m),i,m,a(jt),b(jt))
+ tng(jt)=a(jt)*b(jt).gt.0
+ end do
+ do m=i+1,k
+ if(m.eq.j)cycle
+ if(m.eq.l)cycle
+ jt=(m-1)*(m-2)/2+i
+ call findab(t_in(1:n,m),t_in(1:n,i),m,i,a(jt),b(jt))
+ tng(jt)=a(jt)*b(jt).gt.0
+ end do
+ end subroutine update_pairs3
+
+ subroutine hexagon
+ logical(kind=1),dimension(3,3):: c
+ integer(kind=4),dimension(3) :: v
+ integer(kind=4) :: i,j,m
+ nhex=0
+ do i=1,n-2
+ v(1)=i
+ do j=i+1,n-1
+ if(t_in(i,j).eqv.t_in(j,i))cycle
+ v(2)=j
+ do m=j+1,n
+ v(3)=m
+ c=t_in(v,v)
+ if(all(c.eqv.hexa).or.all(c.eqv.hexb))then
+ nhex=nhex+1
+ hexmat(1:3,nhex)=v
+ endif
+ end do
+ end do
+ end do
+ end subroutine hexagon
+
+ subroutine rand1(x)
+ ! see Ripley, B.D., Stochastic Simulation. New-York:Wiley, 1987, pp. 37-39. (generator 4)
+ integer(kind=4) x,p,q,r,a,b,c
+
+ data a,b,c,p/127773,16807,2836,2147483647/ ! generator 4
+ q=x/a
+ r=mod(x,a)
+ x=b*r-c*q
+ if(x.lt.0)x=x+p
+ end subroutine rand1
+
+ subroutine rand_integer42(iu,a,x1,x2,k)
+ ! draw a uniformly distributed integer from {1, 2, ..., A}
+ ! the rejection rate is (P - BOUND)/P
+ ! where P equals 2**31-2 (= the number of different values RAND1 can take)
+ ! and BOUND equals A*(P/A)
+ ! example: for A = 1000, the rejection rate is 3.008E-7
+ ! Notice that zero as result of the draw leads to rejection (see !***)
+ ! The routine calls RAND2 with the K-th set of coefficients
+
+ integer(kind=4) iu,a,x1,x2,bound,k
+ integer(kind=4),parameter :: p = 2147483646
+ ! P is the number of values X can take: 2**31 - 2,
+ ! because 0 is excluded
+ bound=(p/a)*a
+ do
+ call rand2(x1,x2,k)
+ if(x2.eq.0)cycle
+ if(x2.le.bound)exit !*** LE not LT
+ end do
+ iu=1+mod(x2,a)
+
+ end subroutine rand_integer42
+
+ subroutine rand2(x1,x2,k)
+
+ ! Random number generators MRG (multiple recursive generator)
+ ! Lih-Yuan Deng and Dennis K.J. Lin (2000).Random Number Generation for the New Century,
+ ! The American Statistician, vol 54, no. 2, pp. 145-150
+ ! To compute the formulae, the method in the referenced article is used. B(K) is the table as published
+ ! A(K)=P/B(K) and C(K)=P-A(K)*B(K), P = 2**31 -1
+
+ integer (kind=4)::k
+ integer (kind=4)::x1,x2,p,q,r,y
+ integer (kind=4), dimension(25)::b,c,a
+
+ data b/26403,33236,36673,40851,43693,27149,33986,36848,40961, &
+ 44314,29812,34601,37097,42174,44530,30229,36098,37877, &
+ 42457,45670,31332,36181,39613,43199,46338/
+ data a/81334,64613,58557,52568,49149,79099,63187,58279,52427, &
+ 48460,72034,62064,57888,50919,48225,71040,59490,56696, &
+ 50580,47021,68539,59353,54211,49711,46343/
+ data c/22045, 5979,22786,28279,16390,24896,10265,19055,21300, &
+ 27207, 6039, 7183,12511,25741,24397,15487,13627, 9255, &
+ 8587,34577,19699,32754,23304,18158,41713/
+ data p/2147483647/
+
+ q=x1/a(k)
+ r=mod(x1,a(k))
+ y=b(k)*r-c(k)*q
+ if(y.ge.x2-p)then
+ y=y-x2
+ else
+ y=y+(p-x2)
+ endif
+ if(y.lt.0)y=y+p
+ x1=x2
+ x2=y
+
+ end subroutine rand2
+
+ end subroutine sampler
+
+
+
+ subroutine unpack(vec,words_per_row,t_out,n,k)
+
+
+ integer(kind=4) offset,words_per_row,n,k
+ integer(kind=4),dimension(n*words_per_row)::vec
+ integer(kind=4) i,j,it,ib,ie,ioff
+ ! matrix t_in is not needed
+ !logical(kind=1),dimension(n,k) :: t_in
+ integer(kind=4),dimension(n,k) :: t_out
+
+ t_out=0 ! intialize t_out
+ ioff=0
+ do i=1,n
+ ib=1
+
+ do iw=1,words_per_row
+ ioff=ioff+1
+ ie=min(ib+31,k)
+ it=-1
+ do j=ib,ie
+ it=it+1
+ !t_in(i,j)=btest(vec(ioff),it)
+ !replace the preceding statement by
+ if(btest(vec(ioff),it)) t_out(i,j)=1
+ end do
+ ib=ie+1
+ end do
+
+ end do
+
+ end subroutine unpack
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-raschsampler.git
More information about the debian-science-commits
mailing list