"Erdos-Renyi", Parameters => (n,D,p), Generate => ()->randomMonomialSet(R,D,p) } ) ER (PolynomialRing,ZZ,RR) := (R,D,p) -> ( if p<0.0 or 1.0

"Erdos-Renyi",
Parameters => (R,D,p),
Generate => ()->randomMonomialSet(R,D,p)
}
)
ER (ZZ,ZZ,ZZ) := (n,D,M) -> (
if n<1 then error "n expected to be a positive integer";
x := symbol x;
R := QQ[x_1..x_n];
new Model from {
Name => "Erdos-Renyi",
Parameters => (n,D,M),
Generate => ()->randomMonomialSet(R,D,M)
}
)
ER (PolynomialRing,ZZ,ZZ) := (R,D,M) -> (
new Model from {
Name => "Erdos-Renyi",
Parameters => (R,D,M),
Generate => ()->randomMonomialSet(R,D,M)
}
)
ER (ZZ,ZZ,List) := (n,D,pOrM) -> (
if n<1 then error "n expected to be a positive integer";
if #pOrM != D then error "pOrM expected to be a list of length D";
if not all(pOrM, q->instance(q, ZZ)) and not all(pOrM, q->instance(q,RR))
then error "pOrM must be a list of all integers or all real numbers";
x := symbol x;
R := QQ[x_1..x_n];
new Model from {
Name => "Erdos-Renyi",
Parameters => (n,D,pOrM),
Generate => ()->randomMonomialSet(R,D,pOrM)
}
)
ER (PolynomialRing,ZZ,List) := (R,D,pOrM) -> (
if #pOrM != D then error "pOrM expected to be a list of length D";
if not all(pOrM, q->instance(q, ZZ)) and not all(pOrM, q->instance(q,RR))
then error "pOrM must be a list of all integers or all real numbers";
if all(pOrM, q->instance(q,RR)) and any(pOrM,q-> q<0.0 or 1.0 "Erdos-Renyi",
Parameters => (R,D,pOrM),
Generate => ()->randomMonomialSet(R,D,pOrM)
}
)
sample = method(TypicalValue => Sample)
sample (Model, ZZ) := (M,N) -> (
if N<1 then stderr << "warning: N expected to be a positive integer" << endl;
s:=new Sample;
s.ModelName = M.Name;
s.Parameters = M.Parameters;
s.SampleSize = N;
s.Data = apply(N,i->M.Generate());
s
)
sample String := filename -> (
if not isDirectory filename then error "expected a directory";
modelFile := realpath filename | "Model.txt";
model := lines read openIn modelFile;
s := new Sample;
s.ModelName = model#1;
s.Parameters = value toString stack drop(model,{0,1});
s.SampleSize = value model#0;
dataFile := realpath filename | "Data.txt";
s.Data = value read openIn dataFile;
s
)
getData = method()
getData Sample := s -> (s.Data)
writeSample = method()
writeSample (Sample, String) := (s, dirname) -> (
if fileExists dirname then (
stderr << "warning: directory or file with this name already exists." << endl;
if not isDirectory dirname then (
stderr << "warning: overwrting file." << endl;
removeFile dirname;
mkdir dirname;
);
if isDirectory dirname then (
stderr << "warning: updating name of directory:" << endl;
dirname = dirname |"New";
stderr << dirname << endl;
mkdir dirname;
);
)
else mkdir dirname;
realpath dirname | "Model.txt" << s.SampleSize << endl << s.ModelName << endl << serialize s.Parameters << close;
realpath dirname | "Data.txt" << serialize s.Data << close; -- Write other data
)
statistics = method(TypicalValue => HashTable)
statistics (Sample, Function) := HashTable => (s,f) -> (
fData := apply(getData s,f);
histogram := tally fData;
if not(class fData_0 === ZZ) then (
-- in case the data is actually a BettiTally type, then we are able to get the mean and stddev of the tables:
if (class fData_0 === BettiTally) then (
-- compute the average (entry-wise) tally table:
dataSum := sum histogram;
dataMean := mat2betti(1/s.SampleSize *(sub(matrix(dataSum), RR)));
-- compute the standard deviation (entry-wise) of the Betti tables:
dataMeanMtx := matrix dataMean;
dataVariance := 1/s.SampleSize * sum apply(fData, currentTally -> (
mtemp := new MutableMatrix from dataMeanMtx;
currentTallyMatrix := matrix currentTally;
apply(numrows currentTallyMatrix, i->
apply(numcols currentTallyMatrix, j->
(
--compute mtemp_(i,j) := (bMean_(i,j) - bCurrent_(i,j)):
mtemp_(i,j) = mtemp_(i,j) - currentTallyMatrix_j_i
)
)
);
--square entries of mtemp, to get (bMean_(i,j) - bCurrent_(i,j))^2:
mtemp = matrix pack(apply( flatten entries mtemp,i->i^2), numcols mtemp)
)
);
-- dataStdDev := dataVariance^(1/2); -- <--need to compute entry-wise for the matrix(BettyTally)
dataStdDev := mat2betti matrix pack(apply( flatten entries dataVariance,i->sqrt i), numcols dataVariance);
new HashTable from {
Mean=>mat2betti dataMeanMtx,
StdDev=>dataStdDev,
Histogram=>histogram
}
) else (
stderr << "Warning: the statistics method is returning only the Tally of the outputs of
your function applied to the sample data. If you want more information, such as mean and
standard deviation, then ensure you use a function with numerical (ZZ) or BettiTally output." <

o -> (R,D,p,N) -> ( if p<0.0 or 1.0

o -> (n,D,M,N) -> ( if N<1 then stderr << "warning: N expected to be a positive integer" << endl; if not (instance(o.VariableName,Symbol) or instance(o.VariableName,String) or instance(o.VariableName,IndexedVariableTable)) then error "expected VariableName to be a string or symbol"; x := toSymbol o.VariableName; R := o.Coefficients[x_1..x_n]; apply(N,i-> randomMonomialSet(R,D,M,o)) ) randomMonomialSets (PolynomialRing,ZZ,ZZ,ZZ) := List => o -> (R,D,M,N) -> ( if N<1 then stderr << "warning: N expected to be a positive integer" << endl; apply(N,i-> randomMonomialSet(R,D,M,o)) ) randomMonomialSets (ZZ,ZZ,List,ZZ) := List => o -> (n,D,pOrM,N) -> ( if n<1 then error "n expected to be a positive integer"; if N<1 then stderr << "warning: N expected to be a positive integer" << endl; if not (instance(o.VariableName,Symbol) or instance(o.VariableName,String) or instance(o.VariableName,IndexedVariableTable)) then error "expected VariableName to be a string or symbol"; x := toSymbol o.VariableName; R := o.Coefficients[x_1..x_n]; apply(N,i-> randomMonomialSet(R,D,pOrM,o)) ) randomMonomialSets (PolynomialRing,ZZ,List,ZZ) := List => o -> (R,D,pOrM,N) -> ( if N<1 then stderr << "warning: N expected to be a positive integer" << endl; apply(N,i-> randomMonomialSet(R,D,pOrM,o)) ) randomMonomialSet = method(TypicalValue => List, Options => { Coefficients => QQ, VariableName => "x", Strategy => "ER" }) randomMonomialSet (ZZ,ZZ,RR) := List => o -> (n,D,p) -> ( if p<0.0 or 1.0

o -> (R,D,p) -> ( if p<0.0 or 1.0

o -> (n,D,M) -> (
if n<1 then error "n expected to be a positive integer";
if not (instance(o.VariableName,Symbol) or instance(o.VariableName,String) or instance(o.VariableName,IndexedVariableTable)) then error "expected VariableName to be a string or symbol";
x := toSymbol o.VariableName;
R := o.Coefficients[x_1..x_n];
randomMonomialSet(R,D,M)
)
randomMonomialSet (PolynomialRing,ZZ,ZZ) := List => o -> (R,D,M) -> (
if M<0 then stderr << "warning: M expected to be a nonnegative integer" << endl;
if o.Strategy === "Minimal" then error "Minimal not implemented for fixed size ER model";
allMonomials := flatten flatten apply(toList(1..D),d->entries basis(d,R));
C := take(random(allMonomials), M);
if C==={} then {0_R} else C
)
randomMonomialSet (ZZ,ZZ,List) := List => o -> (n,D,pOrM) -> (
if n<1 then error "n expected to be a positive integer";
if not (instance(o.VariableName,Symbol) or instance(o.VariableName,String) or instance(o.VariableName,IndexedVariableTable)) then error "expected VariableName to be a string or symbol";
x := toSymbol o.VariableName;
R := o.Coefficients[x_1..x_n];
randomMonomialSet(R,D,pOrM,o)
)
randomMonomialSet (PolynomialRing,ZZ,List) := List => o -> (R,D,pOrM) -> (
if #pOrM != D then error "pOrM expected to be a list of length D";
if not all(pOrM, q->instance(q, ZZ)) and not all(pOrM, q->instance(q,RR))
then error "pOrM must be a list of all integers or all real numbers";
B := {};
if all(pOrM,q->instance(q,ZZ)) then (
if o.Strategy === "Minimal" then (
currentRingM := R;
apply(D, d->(
chosen := take(random(flatten entries basis(d+1, currentRingM)), pOrM_d);
B = flatten append(B, chosen/(i->sub(i, R)));
currentRingM = currentRingM/promote(ideal(chosen), currentRingM)
)
)
) else B = flatten apply(toList(1..D), d->take(random(flatten entries basis(d,R)), pOrM_(d-1)));
) else if all(pOrM,q->instance(q,RR)) then (
if any(pOrM,q-> q<0.0 or 1.0(
chosen := select(flatten entries basis(d+1, currentRing), m->random(0.0,1.0)<=pOrM_d);
B = flatten append(B, chosen/(i->sub(i, R)));
currentRing = currentRing/promote(ideal(chosen), currentRing)
)
)
) else B = flatten apply(toList(1..D),d-> select(flatten entries basis(d,R),m-> random(0.0,1.0)<=pOrM_(d-1)));
);
B = apply(B,m->sub(m,R));
if B==={} then {0_R} else B
)
bettiStats = method(TypicalValue =>Sequence, Options =>{IncludeZeroIdeals=>true, SaveBettis => "", CountPure => false, Verbose => false})
bettiStats List := o-> (ideals) -> (
N := #ideals; Z:=0;
if o.SaveBettis != "" then (
if fileExists o.SaveBettis then (
stderr << "warning: filename already exists. Overwriting." << endl;
removeFile o.SaveBettis;
);
);
if not o.IncludeZeroIdeals then (
(ideals,Z) = extractNonzeroIdeals(ideals);
if o.Verbose then stdio << "There are "<