use application 'ideal';


# Create a random (possible) Wronski System
# Usage: random_wronski(M,lambda,{seed=>r, triangulation=>T, param=>s, ring=>R});
# @param Matrix<Int> M points of the Triangulation (Homogenized)
# @param Vector<Int> lambda lifting function
# @option Int seed for randomizer
# @option topaz::GeometricSimplicialComplex triangulation the Triangulation
# @option Scalar<Rational> param the additional parameter for creating a wronski system
# @option Ring ring The Ring
sub random_wronski($$;$)
{
    # reading parameters
    my ($M, $lambda, $options) = @_;
    my ($n, $d) = ($M->rows(), $M->cols()-1);
    
    # randomizer options
    my $offset = -50; my $range = 100;
    my $seed = $options->{"seed"};
    if (defined($seed)){ srand($seed);};

    # triangulation
    my $T = $options->{"triangulation"};
    if (!defined($options->{"triangulation"})) {
      $T=new topaz::SimplicialComplex(FACETS=>polytope::regular_subdivision($M,$lambda));
    }

    # ring
    my $R = $options->{"ring"};

    # radomized coefficients
    my $coeffs = new Array<Array<Rational> >($d);
    for (my $i=0; $i<$d; ++$i){
	my $temp = new Array<Rational>($d+1);
	map { $temp->[$_] = int(rand($range)) + $offset } (0..$d);
	$coeffs->[$i] = $temp;
    }

    # the additional parameter
    my $s = $options->{"param"};
    if (!defined($options->{"param"})){
	$s=int(rand(99)+1)/100;
    }

    # creating the system
    return wronski_system($M,$lambda,$coeffs,$s,{triangulation=>$T, ring=>$R});
}


# Solves many random wronski systems and checks how many real solutions there are.
# the data will be stored in an output file
# All systems are seeded, so that one can check the individual systems
# Usage: solve_random_wronski(M,lambda,n,{seedoffset=>r, param=>s, upperparam=>s0, file=>"path", triangulation=>T});
# @param Matrix<Int> M points of the Triangulation (Homogenized)
# @param Vector<Int> lambda lifting function of the triangulation
# @param Int n Number of Systems to solve
# @option Int seedoffset the starting seed
# @option Scalar<Rational> param the additional parameter of the wronski system
# @option Scalar<Rational> upperparam the upper bound for the additional parameter
# @option File the output file of the analysis
# @option topaz::GeometricSimplicialComplex T Triangulation
sub solve_random_wronski($$$;$)
{
    # reading the parameters
    my ($M, $lambda, $n, $options) = @_;

    # triangulation
    my $T = $options->{"triangulation"};
    if (!defined($options->{"triangulation"})) {
	$T=new topaz::SimplicialComplex(FACETS=>polytope::regular_subdivision($M,$lambda));
    }

    # output file
    my $file = $options->{"file"};
    if (!defined($options->{"file"})){
	$file="solve_random_wronski.csv";
    }

    # upper bound for the additional parameter
    my $s0 = $options->{"upperparam"};
    if (!defined($options->{"upperparam"})){
	$s0 = 1;
    }

    # seed start
    my $seedoffset = $options->{"seedoffset"};
    if (!defined($options->{"seedoffset"})){
	$seedoffset=0;
    }

    # Notification
    print "\n====================\n";
    print "==== Test Start ====\n";
    print "====================\n";

    # Kushnirenko's bound  
    my $p = new polytope::Polytope(POINTS=>$M);
    my $max_sol = fac($p->DIM)*$p->VOLUME;

    # Number of systems which are not Wronski Systems
    my $no_wronski = 0;
    
    # Open output file
    open OUTPUT_FILE, ">$file"
    or die "can't create outputfile $file: $!";
    print OUTPUT_FILE '"System";"Seed";"Parameter s=";"Number of real roots"'. "\n";
    $|=1;

    # Iterating through the systems
    for (my $seed=$seedoffset + 1; $seed<=$seedoffset+$n; ++$seed){
	# the additional parameter
	my $s = $options->{"param"};
	if (!defined($options->{"param"})){
	    $s=int(rand($s0*100-1)+1)/100;
	}

	# creating a (possible) wronski system (with seed according to the iteration number)
	my $randws = random_wronski($M, $lambda,{param=>$s,seed=>$seed,triangulation=>$T});

	# Solving the system
	my $roots;
	eval{$roots = $randws->SOLVE;};

	# if the system is not 0-dimensional then it is not a wronski system
	# and we forget that we solved it.
	if($@){
	    print "(".($seed-$no_wronski).") Seed: ".$seed." -- not 0-dim!\n";
	    ++$no_wronski;
	    ++$n;
	    next;
	}

	# if Kushnirenko's bound isn't met it is not a wronski system
	# note: this is the case if some crucial coefficiants are 0
	# note2: we check for "<" instead of "=" because 0 could be a solution
	if ($roots->rows() < $max_sol){
	    print "(".($seed-$no_wronski).") Seed: ".$seed." -- no wronski!  (expected: $max_sol   got: ". $roots->rows().")\n";
	    ++$no_wronski;
	    ++$n;
	    next;
	}

	# counting the real solutions
	my $real_sols = 0;
	for (my $i = 0; $i<$roots->rows(); ++$i){
	    my $is_real = 1;
	    for (my $j = 0; $j<$roots->cols(); ++$j){
		if ($roots->[$i]->[$j]->[1] != 0) {
		    $is_real = 0; 
		    last;
		}
	    }
	    if ($is_real == 1) {++$real_sols;};
	}

	# writing everything to the outputfile
	print OUTPUT_FILE '"'.($seed-$no_wronski).'";"'.$seed.'";"'.$s.'";"'.$real_sols.'"'."\n";


	# notification
	if ($seed % 10000 == 0) {
	    print "($seed) Counter: ".($seed - $no_wronski)." -- No Wronski: $no_wronski\n";
	}
    }

    close OUTPUT_FILE;
    print "==== Test End ====\n";
}