(* :Title: JRandomArray *)

(* :Author: Mark Fisher *)

(* :Version: 1.0 December 2004 *)

(* :Mathematica Version: 5.0 *)

(* :Requirements:
This package requires colt.jar and JRandomArray.class.
*)

(* :Summary:
JRandomArray uses JLink to implement a variety of pseudo random number
generators from the Java colt distribution
(see http://hoschek.home.cern.ch/hoschek/colt/index.htm).

By computing arrays of random numbers directly in (a small amount of) Java
code itself, the overall speed is competitive with Mathematica itself.
(See RandomArray.java for the requisite Java.)

One should be aware that with the exception of MersenneTwister64, the
random number engines compute only 32 bits. (The remaining bits are filled
in as zeros.) MersenneTwister64 returns about 54 bits (equivalent to about
16 decimals). I suspect this truncation is produced by the Java
implementation of the random engine running under Windows: Returning the
random doubles as BigDecimals leaves the number of significant bits
unchanged. By the way, Mathematica's Random[] does not pad with zeros
beyond 54 bits but instead returns some sort of noise out to about 62 bits.
*)

SeedRandomEngine::usage = "SeedRandomEngine[randomEngineClass] returns a
random engine object from the specified class using
JavaNew[\"java.util.Date\"] as the default seed. Alternatively, one can
specify the seed explicitly using JavaNew[randomEngineClass, seed], where
seed is a positive integer or Java date object. See the usage notes for
RandomEngineClass for a list of valid random engine class names. Example:\n\n
engineObj = SeedRandomEngineClass[\"MersenneTwister64\"]."

InitializeDistribution::usage = "InitializeDistribution[distributionClass,
randomEngineObject][pars] returns a distribution object from the specified
class driven by the given random engine object and using the specified
parameters. For example,\n\n
normalObj = InitializeDistribution[\"Normal\", engineObj][0,1],\n\n
where mersenneTwister64 is a seeded random engine object.
To reseed a distribution object, one must first reseed the
random engine and then reinitialize the distribution. Each distribution
has a method to change the parameter values. For example, to change the
standard deviation of the normal distribution from 1 to 3, use\n\n
normalObj@setState[0, 3].\n\n
See the usage notes for DistributionClass for a list of valid distribution
class names and parameter reset methods."

JRandomArray::usage = "JRandomArray[distributionObject, n] returns a vector
of n pseudo random numbers generated according to the distribution object.
JRandomArray[distributionObject, {n1, n2, ...}] returns a higher
dimensional array of random numbers. JRandomArray[distributionObject]
returns a single element."

(* setup JLink *)
Needs["JLink`"]
InstallJava[]

(* 	To generate 10^7 random variables, increase available memory:
	to allocate 50MB to start (not necessary) and allocate 200MB.
	This has the side effect (on Windows) of opening a cmd window
	which remains until UninstallJava[] is executed or the kernel
	is quit.

InstallJava[CommandLine -> "java -Xms50m -Xmx200m"]
*)


(* make the classes accessable: colt.jar and RandomArray.class *)
(**** you should modify this path to suit your directory structure ****)

AddToClassPath["d:/data/packages/java"]

(* JRandomArray class implements the static method jRandomArray
	that fills a vector with random numbers and returns it *)
randomArrayClass = LoadJavaClass["JRandomArray"]

SeedRandomEngine[randomEngineClass_String] :=
	JavaBlock[
	JavaNew[
		RandomEngineClass[randomEngineClass],
		JavaNew["java.util.Date"]]
	]

InitializeDistribution[distributionClass_String,
		randomEngineObj_?JavaObjectQ][pars__?NumberQ] :=
	JavaNew[
		DistributionClass[distributionClass],
		pars,
		randomEngineObj]

JRandomArray[distributionObject_?JavaObjectQ, n_Integer?Positive] :=
	JRandomArray`jRandomArray[distributionObject, n]

JRandomArray[distributionObject_?JavaObjectQ, n:{__Integer?Positive}] :=
	With[{array = JRandomArray[distributionObject, Times @@ n]},
	If[Length[n] == 1,
		(* then *)
		Partition[array, 1],
		(* else *)
		Fold[Partition, array, Drop[Reverse[n], -1]]
		]
	]

JRandomArray[distributionObject_?JavaObjectQ] :=
	JRandomArray[distributionObject, 1][[1]]

(* helpful lists *)
shortRandomEngineClassNames = {
	"DRand", "MersenneTwister", "MersenneTwister64",
	"Ranecu", "Ranlux", "Ranmar"}

MapThread[Set, {
	RandomEngineClass /@ shortRandomEngineClassNames,
	LoadJavaClass /@ {
		"cern.jet.random.engine.DRand",
		"cern.jet.random.engine.MersenneTwister",
		"cern.jet.random.engine.MersenneTwister64",
		"edu.cornell.lassp.houle.RngPack.Ranecu",
		"edu.cornell.lassp.houle.RngPack.Ranlux",
		"edu.cornell.lassp.houle.RngPack.Ranmar"
	}
	}]

(* need to provide documentation regarding valid parameter values *)
shortDistributionClassNames = {
	"Beta", "Binomial", "ChiSquare", "Exponential",
	"Gamma", "HyperGeometric", "Logarithmic", "NegativeBinomial",
	"Normal", "Poisson", "StudentT", "Uniform"}

MapThread[Set, {
	DistributionClass /@ shortDistributionClassNames,
	LoadJavaClass["cern.jet.random." <> #]& /@ shortDistributionClassNames
	}]

parstable = ToString @
	TableForm[Transpose[{
		shortDistributionClassNames,
		StringDrop[Select[Methods[DistributionClass[#], Inherited -> False][[1]],
			StringMatchQ[#, "*set*"] &][[1]], 5] & /@ shortDistributionClassNames}],
	TableSpacing -> {0, 3}]

DistributionClass::usage = "DistributionClass[className] returns
the JLink object that represents the named distribution class.
Valid distribution classes are the strings given
in first column of the following list:\n" <> parstable <> "\nThe second
column provides the name of the method used to reset the parameter values
and also gives the number and type of parameters for each class."

RandomEngineClass::usage = "RandomEngineClass[className] returns
the JLink object that represents the named random engine class.
Valid class names are the strings\n" <>
	ToString[TableForm[shortRandomEngineClassNames, TableSpacing -> {0}]]

(* :Example: *)
(* use 64 bit twister engine; returns only about 54 significant bits *)
engineObj = SeedRandomEngine["MersenneTwister64"]

(* for uniform (0,1) and normal (0,1) distributions *)
uniformObj = InitializeDistribution["Uniform", engineObj][0, 1]
normalObj = InitializeDistribution["Normal", engineObj][0, 1]

Null; (* don't return anything *)