Skip to content

Commit

Permalink
draft for 2x2 example done
Browse files Browse the repository at this point in the history
  • Loading branch information
rfabbri committed Jun 15, 2024
1 parent 974ba61 commit 1327f89
Show file tree
Hide file tree
Showing 3 changed files with 140 additions and 20 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -180,3 +180,4 @@ ys3d.tmp
ysnew.tmp
scripts/synthdata/results-synth/
tests/Fastor/
original
135 changes: 127 additions & 8 deletions tutorial/end-2x2.m2
Original file line number Diff line number Diff line change
Expand Up @@ -39,15 +39,134 @@
--

restart
needs "parser.m2"
elapsedTime needs "chicago.m2"
needs "MinusUtility.m2"

(p0,sols) = readStartSys "startSys";
-- Pro -------------------------------------------------------------------------
--
-- Set online tracker options here
--
-- null indicates default value
-- TODO: show how to get default values here
-- TODO: are these the same for start and end sys?
-- scan({CorrectorTolerance=>1e-4,
-- EndZoneFactor=>2e-1,
-- InfinityThreshold => 1e6,
-- maxCorrSteps => 3,
-- noOutput => true,
-- numberSuccessesBeforeIncrease => 2,
-- Precision => null,
-- Predictor => RungeKutta4,
-- stepIncreaseFactor => 2,
-- tStep => 5e-2,
-- tStepMin => 1e-5
-- },
-- opt -> setDefault(opt)) -- setDefault is in MonodromySolver
--
-- setDefault(CorrectorTolerance=>1e-8)


-- Noob and Pro ----------------------------------------------------------------
-- Read start system
(p0,sols0) = readStartSys "startSys";

-- Pro -------------------------------------------------------------------------
-- verify options are set as desired
netList((keys options trackHomotopy)/(opt ->(opt, getDefault opt)))
NAGtrace 3
setRandomSeed 0
-- netList((keys options trackHomotopy)/(opt ->(opt, getDefault opt)))
-- NAGtrace 3
-- setRandomSeed 0
--
-- Reads from file
-- p1 = data2parameters("target_system_data.txt');

p0 := if (numcols p0 == 1) then p0 else transpose p0;

-- Noob ------------------------------------------------------------------------
p1 = {2, 3, 4, 6, -2, 11.4};
x1gt = {} -- ground truth solutions of target system

-- p1 = {a, b, c, d, e, f};
-- where {a*(x^2+y^2)+b*x+c, d*x+e*y+f}
--
-- Pro -------------------------------------------------------------------------
-- Convert your physical data to the parameters here.
-- p1 := data2parameters targetData; -- XXX just put any abcdef here

P01 = p0||p1;

-- Pro -------------------------------------------------------------------------
-- P01 := (gammify p0)||(gammify p1); -- include Noob and Pro gammify

-- (pLines, x) = parseFile(1,Tests=>true);
-- elapsedTime K = solveChicago(p0, pLines, sols);
-- Noob and Pro ----------------------------------------------------------------
H01 := specialize(PH, P01);
trackHomotopy(H01, sols0, K)

-- Evaluate check
--
F = xxx
evaluate(F,x0||p0)
evaluate(F,x1||p1)

-- Pro -------------------------------------------------------------------------
-- J0 = evaluate(J,x0||p0); -- Evaluates Jacobian if desired

-- Pro -------------------------------------------------------------------------
-- Example of technique to improve speed and robustness
--
-- -*
-- (p,x)=fabricateChicago CC
-- evaluate(F,x||p)
--
-- (p0,x0) = fabricateChicago(CC)
-- evaluate(F,x0||p0)
-- --J0 = evaluate(J,x0||p0);
-- --272, 279, 286
-- S = subsets(4,2)
-- for s in S do (
-- << s << endl;
-- Jpivots = flatten for i from 0 to 4 list {4*i+s#0,4*i+s#1};
-- F' = F^(Jpivots);
-- J'=diff(matrix{cameraVars},F');
-- J'0 = evaluate(J', x0||p0);
-- << J'0 << endl;
-- S=first SVD J'0;
-- << toString((max S)/(min S)) << endl;
-- << S << endl;
-- )
--
-- norm evaluate(F',x0||p0)
-- transpose x0^{6..11} * p0^{39..44}
-- p0_(45,0)
--
-- netList rsort for k from 272 to 272+15-1 list (
-- --try 289?
-- Jpivots = {2, 3, 6, 7, 10, 11, 14, 15, 18, 19,k,numrows F-3,numrows F-2,numrows F-1};--rowSelector J0;
-- F' = F^(Jpivots);
-- J'=diff(matrix{cameraVars},F');
-- S = first SVD evaluate(J', x0||p0);
-- c = (max S)/(min S);
-- e = evaluate(diff(matrix{cameraVars},F^{k}),x||p);
-- n = norm e;
-- (c, k, n, e)
-- )
-- *-
-- -- INDICES FOR SQUARE SUBSYTEM
-- Jpivots = flatten(for i from 0 to 4 list {4*i+CLBlocks#0,4*i+CLBlocks#1}) | {274,numrows F-3,numrows F-2,numrows F-1};
--
-- assert(# Jpivots==nvars)
-- F' = F^(Jpivots);
-- J'=diff(matrix{cameraVars},F');
-- elapsedTime GS = gateSystem(gateMatrix{dataParams},gateMatrix{cameraVars},F');
-- elapsedTime PH = parametricSegmentHomotopy GS;
-- setDefault(CorrectorTolerance=>1e-8)
--
-- -- filter path jumps during monodromy
-- filterEval = (p,x) -> (
-- -- false iff residual small
-- resid := norm evaluate(F,x||p);
-- -- << "residual: " << resid << endl;
-- (resid > 1e-4)
-- )
--
-- seeds = {16314, 62236, 77740, 69544, 11665, 1055}
-- H = hashTable apply(subsets(4,2),seeds, (ss,s) -> ss=>s)
--
24 changes: 12 additions & 12 deletions tutorial/start-2x2.m2
Original file line number Diff line number Diff line change
Expand Up @@ -37,15 +37,15 @@
restart
needsPackage "SLPexpressions"
needsPackage "MonodromySolver"
needs "MinusUtility.m2" -- put cCode extesion here
-- "gateSystem" exists only in M2 v 1.14
needs "MinusUtility.m2"

-- Noob---------------------------------------------------------------------------
-- Pro doesn't use declareVariable
variables = declareVariable \ {x,y}
params = declareVariable \ {a,b,c,d,e,f}

-- Noob and Pro ------------------------------------------------------------------
-- gateSystem exists only in M2 v 1.14
GS = gateSystem(
matrix{params},
matrix{variables},
Expand All @@ -56,19 +56,18 @@ GS = gateSystem(
)
-------------------------------------------------------------------------------


-- Pro optional
-- Seed random seeds for reproducible
-- Pro
-- random seeds for reproducible runs
-- setRandomSeed H#CLBlocks


-- Noob and Pro ------------------------------------------------------------------
cameraVars = flatten entries vars GS
symbols = flatten entries vars GS
PH = parametricSegmentHomotopy GS

-- Pro only --------------------------------------------------------------------
-- Pro ---------------------------------------------------------------------------

-- YOU set tracker options here
-- Set Monodromy tracker options here
--
-- null indicates default value
-- TODO: show how to get default values here
scan({CorrectorTolerance=>1e-4,
Expand All @@ -87,7 +86,7 @@ scan({CorrectorTolerance=>1e-4,

setDefault(CorrectorTolerance=>1e-8)

-- Pro Gammify-style Randomization ---------------------------------------------------
-- Pro Gammify-style Randomization -------------------------------------------------


-- YOU: adapt for your problem
Expand Down Expand Up @@ -115,11 +114,11 @@ setDefault(CorrectorTolerance=>1e-8)
-- HxHt
h=cCode(
transpose(PH.GateHomotopy#"Hx"|PH.GateHomotopy#"Ht"),
gateMatrix{cameraVars|{PH.GateHomotopy#"T"}|flatten entries PH#Parameters}
gateMatrix{symbols|{PH.GateHomotopy#"T"}|flatten entries PH#Parameters}
)

-- HxH
h=cCode(transpose(PH.GateHomotopy#"Hx"|PH.GateHomotopy#"H"),gateMatrix{cameraVars|{PH.GateHomotopy#"T"}|flatten entries PH#Parameters})
h=cCode(transpose(PH.GateHomotopy#"Hx"|PH.GateHomotopy#"H"),gateMatrix{symbol|{PH.GateHomotopy#"T"}|flatten entries PH#Parameters})

-- Maybe useful
-- cCode PH
Expand All @@ -128,6 +127,7 @@ h=cCode(transpose(PH.GateHomotopy#"Hx"|PH.GateHomotopy#"H"),gateMatrix{cameraVar
-- monodromy needs an initial pair of parameter, solution
-- the command below won't work for most use cases
-- ie) need
-- XXX
(p0, x0) = createSeedPair GS

-- Pro 1 -----------------------------------------------------
Expand Down

0 comments on commit 1327f89

Please sign in to comment.