restart

{*Macaulay2: software system devoted to supporting 
research   in algebraic geometry and commutative algebra*}
--The Authors: Dan Grayson and Mike Stillman
--Goal: Describe algebraic varieties. 


---PART 1-----M2 basics and help
R=QQ[x,y,z]--This is a ring
I=ideal(x^2+y^2-1,x-z)--This is an ideal generated 
degree I 
dim I
codim I
--M2 has very good documentation. 
help degree
help ideal
help apropos
apropos "ide"


---PART 2-----Twisted Cubic
--Elimination
R=QQ[x0,x1,x2,x3]**QQ[s,t]--This ring organizes the variables into two groups. 
degree (x0-s^2)
I1=ideal({x0,x1,x2,x3}-{s^3,s^2*t,s*t^2,t^3})
I2=eliminate({s,t},I1)

--Minors of matrix
R=QQ[x0,x1,x2,x3]
M=matrix{{x0,x1,x2},{x1,x2,x3}}
I=minors(2,M)
help genericMatrix 
--hilbertPolynomial I
--resI=resolution I
class I
theGenerators=I_*---this is a list
class theGenerators
generators I--This is a matrix
flatten entries gens I
J=ideal drop(theGenerators,1)
help drop
numgens I
numgens J
decJ=decompose J
netList decJ--Prints the list in a nice way.
decJ/degree --degree of each element in the list. 
help apply

---PART 3-----Resultants
viewHelp "Macaulay2"
viewHelp "Package"--See the many packages. 
loadPackage"Resultants";
F = genericPolynomials(QQ,{4,4})
first F
time resF = Resultant F
degree resF
--F = genericPolynomials(QQ,{4,4,4})
--time resF = Resultant F
--degree resF

---PART 4-----Other fields
--We can also work over finite fields and multigraded rings. 
kk=ZZ/30103
isPrime 30103
R=kk[x0,x1]**kk[y0,y1]
f= x0^2*y1-x1^2*y0
degree f
multidegree ideal f
--Other fields that are commonly used in M2. 
class QQ
class RR
class (ZZ/30103)
class CC

---------------------------------------------------
--Numerical algebraic geometry
---------------------------------------------------
---PART 5-----NumericalAlgebraicGeometry
kk=CC
R=kk[x1,x2]
--Anton Leykin
--Robert Krone
loadPackage("NumericalAlgebraicGeometry",Reload=>true)
help solveSystem 
theSols=solveSystem({x1^2+x2^2-1,x1-x2})
(sqrt 2)/2
onePoint= first theSols
class onePoint
keys onePoint
onePoint#SolutionStatus
netList pairs onePoint

--These functions allow us to test timing and residuals. 
testTiming=(numVars,theDegree)->(
  R=kk[apply(numVars,i->"x"|i+1)];
--  F:=apply(numVars,i->random(theDegree,R)+random CC);
  F=apply(numVars,i->sum apply(theDegree+1,d->random(d,R)));
  startTime:=currentTime();
  win=solveSystem(F);
  endTime:=currentTime();
  totalTime:=(endTime-startTime);
  timePerPath:=totalTime/#win*1.0;
  ("Total Time: ",totalTime,
      "Time per path", timePerPath,
      "Bezout Bound",theDegree^numVars,
      "Found solutions",#win))
testTiming(2,40)

testResiduals=(k,aSol)->(
    thePoly:=F_k;
    theSubs:=apply(gens ring F_k,aSol#Coordinates,(i,j)->i=>j);
    sub(F_k,theSubs)
    )

testResiduals(1,first win)
theResids=apply(win,i->abs testResiduals(0,i));
max theResids
min theResids

---PART 6-----Bertini.m2
restart
loadPackage"Bertini"
--Jose Israel Rodriguez
--Dan Bates
--Anton Leykin
--Elizabeth Gross 
loadPackage"Resultants"
F = genericPolynomials(QQ,{4,4})
R= ring first F
theCoefficients=gens coefficientRing R
S=CC[s,t] 
coefSubs=apply(theCoefficients,i->i=>random(1,S))
--This makes a Bertini input file on your computer.
makeB'InputFile(storeBM2Files,
    HomVariableGroup=>{gens S,gens R},
    B'Polynomials=>F,
    B'Functions=>coefSubs   
    )
readFile(storeBM2Files,"input",10000)--see the input file
runBertini(storeBM2Files)--This calls Bertini.
readFile(storeBM2Files)--This reads a text file summarizing the results of the homotopy computation.
readFile(storeBM2Files,"nonsingular_solutions",10000)
--This imports just the coordinates.
theWitnessPointsCoordinates=importSolutionsFile(storeBM2Files);
first theWitnessPointsCoordinates
--To get more information regarding condition number and etc, import the main data file.
readFile(storeBM2Files,"main_data",1000)
witnessPoints=importMainDataFile(storeBM2Files);
aWitnessPoint=first witnessPoints
keys aWitnessPoint
netList pairs  aWitnessPoint