ComplexityOneVariety.mm

module ComplexityOneVariety()
    option object;

    # These fields are guaranteed to be filled when a ComplexityOneVariety is created.
    export P;
    export Sigma := undefined;
    export A := undefined;

    # Variables and relations in the Cox Ring
    export variables := undefined;
    export monomials := undefined;
    export relations := undefined;

    ###################################################################################
    ### These fields are only computed when needed. Use the getters below for them. ###
    ###################################################################################

    # The maximal X-cones of the fan Sigma, i.e. the maximal cones of the fan of the minimal
    # ambient toric variety
    local maximalXCones := undefined;

    # The fan of the minimal ambient toric variety, stored as a `FAN` from the convex package
    local minimalAmbientFan := undefined;

    # Says whether the variety is Q-factorial, i.e. if every weil divisor admits a Cartier multiple.
    # This is equivalent to whether the minimal ambient toric variety is Q-factorial.
    local isQfactorialVal := undefined;

    # A list of local picard indices, one for each X-cone. Might contain infinite values
    # if the variety is not Q-factorial.
    local localPicardIndices := undefined;

    # The picard index of the variety, i.e. the index of the picard group in the divisor class group.
    # It is given as the least common multiple of all local picard indices.
    # It is infinity if and only if the variety is not Q-factorial
    local picardIndex := undefined;

    # Says whether the variety is factorial, i.e. if every weil divisor is Cartier
    # This is equivalent to whether the minimal ambient toric variety is smooth.
    # It is also euivalent to the picard index being one.
    local isFactorialVal := undefined;

    # Says whether the variety is Q-gorenstein, i.e. if the anticanonical class admits a Cartier multiple.
    local isQgorensteinVal := undefined;

    # A list of local gorenstein indices, one for each X-cone. Might contain infinite values
    # if the variety is not Q-gorenstein.
    local localGorensteinIndices := undefined;

    # The gorenstein index of the variety.
    # It is given as the least common multiple of all local gorenstein indices.
    # It is infinity if and only if the variety is not Q-gorenstein.
    local gorensteinIndex := undefined;

    # Says whether the variety is gorenstein, i.e. if the anticanonical class is Cartier.
    # This is equivalent to the gorenstein index being one.
    local isGorensteinVal := undefined;

    # The local quotients of picard indices and gorenstein indices
    local localGorensteinQuotients := undefined;

    # The quotient of picard index and gorenstein index
    local gorensteinQuotient := undefined;


    local ampleCone := undefined;

    local isFanoVal := undefined;

    ###############################################################################
    ### The following fields are only defined for K*-surfaces, i.e. when s = 1. ###
    ###############################################################################
    local intersectionMatrix := undefined;
    local anticanonicalSelfIntersection := undefined;

    # This field is only filled if the variety is created as some resolution of singularities from another variety
    export exceptionalDivisorsIndices := undefined;

    export ModuleApply :: static := proc()
        Object(ComplexityOneVariety, _passed);
    end;

    export setCoefficientMatrix :: static := proc(self :: ComplexityOneVariety, A :: Matrix)
        local P, f, lss, i, j, alpha, T, S;
        P := self:-P;
        f := P:-format;
        lss := P:-lss;
        self:-variables := [seq(seq(T[i,j], j = 1 .. f:-ns[i]), i = 0 .. f:-r), seq(S[i], i = 1 .. f:-m)];
        self:-monomials := [seq(mul([seq(T[i,j] ^ lss[i][j], j = 1 .. f:-ns[i])]), i = 0 .. f:-r)];

        alpha := (i,j) -> Determinant(Matrix([Column(A, [i,j])]));
        self:-relations := [seq(alpha(i+1,i+2) * self:-monomials[i] + alpha(i+2, i) * self:-monomials[i+1] + alpha(i, i+1) * self:-monomials[i+2], i = 1 .. f:-r - 1)];
    end proc;

    (*
    This method creates a T-Variety of complexity one from various kinds of data. 
    The standard input method looks like this:

    P :: PMatrix, Sigma :: {set(set(integer)), Vector, list}, A :: Matrix.

    In some cases (see below), `Sigma` is optional. `A` is always optional.

    The parameter `Sigma` encodes the fan of the ambient toric variety. It can either be given directly as a set of
    cones, where each cone is encoded by the set of column indices of P it contains, or it can be given by 
    a weight in the rational vector space associated to the class group of P lying in the moving cone of P. 
    In the latter case, we take the fan to be the gale dual of the bunch of orbit cones generated by that weight.

    The parameter `Sigma` is optional in the following three cases:

    (1) s = 1. In this case, we have a C*-surface and there is only one possible minimal ambient fan, which we
        can explicitly construct from P (see Cox Rings 5.4.1.6).
    (2) classGroupRank = 1. In this case, the number of rays of the fan is one more than the ambient dimension, hence
        there is only one complete fan having P as generator matrix.
    (3) P admits a Fano variety. In this case, the anticanonical class is contained in the moving cone of P, so we can 
        use it to define a bunch of cones and hence an ambient fan. The resulting variety will be the unique Fano variety
        having P as its PMatrix. Note however, that there are other possible non-Fano varieties with P as PMatrix.

    The parameter `A` is the coefficient matrix for the trinomial euqations defining the variety.
    It is an optional parameter. If it is not provided, it will be left undefined. It can always be
    provided later with the `setCoefficientMatrix` procedure.

    *)
    export ModuleCopy :: static := proc(self :: ComplexityOneVariety, proto :: ComplexityOneVariety, P :: PMatrix)
        local numColumns, i, j, k, Sigma, taus, sigma_plus, sigma_minus, taus_plus, taus_minus, w, candidates, minimalBunchCones, cone, A, vplus_index, vminus_index;

        self:-P := P;

        # Check the arguments, set Sigma and A if present.
        if _npassed > 3 then
            if type(_passed[4], Matrix) then
                # No Sigma provided, coefficient Matrix is fourth argument.
                A := _passed[4];
                if RowDimension(A) <> 2 or ColumnDimension(A) <> P:-r + 1 then
                    error "The coefficient matrix must be a (2 x r)-Matrix. Here, r = %1", P:-r;
                end if;
                for i from 1 to P:-r + 1 do
                    for j from i+1 to P:-r + 1 do
                        if Determinant(Matrix([Column(A, [i,j])])) = 0 then
                            error "The columns of the coefficient matrix must be linearly independent. Here the %1-th and the %2-th columns are linearly dependent.", i, j;
                        end if;
                    end do;
                end do;
                self:-A := A;
            elif type(_passed[4], set(set(integer))) then
                # Fan Sigma provided
                if not foldr((c1,c2) -> c1 union c2, {}, op(_passed[4])) = {seq(1 .. P:-numCols)} then
                    error "The cones of `Sigma` must cover {1 .. %1}", P:-numCols;
                end if;
                Sigma := _passed[4];

                self:-Sigma := _passed[4];
            elif type(_passed[4], {Vector, list(integer)}) then
                # Weight w provided. Compute the associated fan.
                w := _passed[4];

                if not containsrelint(getMovingCone(P), w) then
                    error "The given weight w = %1 does not lie in the moving cone of P.", convert(w, list);
                end if;

                # TODO: Make this computation more efficient by incrementally searching only the minimal cones containing w.
                candidates := [op(combinat[powerset]({seq(1 .. P:-n + P:-m)}) minus {{}})];
                minimalBunchCones := {};
                for i from 1 to nops(candidates) do
                    cone := candidates[i];
                    # If we already have a cone that's contained in this one, skip it.
                    if select(c -> c subset cone, minimalBunchCones) = {} then
                        if containsrelint(poshull(Column(getDegreeMatrixFree(P), [op(cone)])), w) then
                            minimalBunchCones := {op(minimalBunchCones), cone};
                        end if;
                    end if;
                end do;
                # We dualize to get the maximal cones of the associated fan.
                self:-Sigma := map(c -> {seq(1 .. P:-n + P:-m)} minus c, minimalBunchCones);
            else
                error "Expected second argument to be of type set(set(integer), Vector, list or Matrix";
            end if;

            if _npassed > 4 then
                if type(_passed[5], Matrix) then
                    A := _passed[5];
                    if RowDimension(A) <> 2 or ColumnDimension(A) <> P:-r + 1 then
                        error "The coefficient matrix must be a (2 x r)-Matrix. Here, r = %1", P:-r;
                    end if;
                    for i from 1 to P:-r + 1 do
                        for j from i+1 to P:-r + 1 do
                            if Determinant(Matrix([Column(A, [i,j])])) = 0 then
                                error "The columns of the coefficient matrix must be linearly independent. Here the %1-th and the %2-th columns are linearly dependent.", i, j;
                            end if;
                        end do;
                    end do;
                    self:-A := A;
                else
                    error "Expected third argument to be of type Matrix";
                end if;
            end if;
        end if;

        if not type(self:-A, undefined) then
            setCoefficientMatrix(self, self:-A);
        end if;

        # If no Sigma has been provided, check if we are in one of the allowed cases (1)-(3)
        # and compute it.
        if type(self:-Sigma, undefined) then
            if P:-s = 1 then
                # Here, X is a C*-surface.
                # In this case, the moving cone equals the semiample cone, hence there is only one
                # possible fan for the P-Matrix. We can write it down explicitly, following Construction 5.4.1.6 of "Cox Rings".
                # Note that we do not require the P-Matrix to be slope-ordered.
                
                sigma_plus := {seq(doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][1]), i = 0 .. P:-r)};
                sigma_minus := {seq(doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][P:-ns[i]]) , i = 0 .. P:-r)};
                
                taus := {seq(seq({
                    doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][j]), 
                    doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][j+1])}, 
                    j = 1 .. P:-ns[i] - 1), i = 0 .. P:-r)};
                
                if P:-case = "EE" then
                    self:-Sigma := {sigma_plus} union taus union {sigma_minus};
                elif P:-case = "PE" then
                    taus_plus := {seq({doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][1]), P:-n + 1} , i = 0 .. P:-r)};
                    self:-Sigma := taus_plus union taus union {sigma_minus};
                elif P:-case = "EP" then
                    taus_minus := {seq({doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][P:-ns[i]]), P:-n + 1} , i = 0 .. P:-r)};
                    self:-Sigma := {sigma_plus} union taus union taus_minus;
                elif P:-case = "PP" then
                    if P:-d[1, P:-n + 1] = 1 then
                        vplus_index := P:-n + 1;
                        vminus_index := P:-n + 2;
                    elif P:-d[1, P:-n + 1] = -1 then
                        vplus_index := P:-n + 2;
                        vminus_index := P:-n + 1;
                    end if;
                    taus_plus := {seq({doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][1]), vplus_index} , i = 0 .. P:-r)};
                    taus_minus := {seq({doubleToSingleIndex(P:-format, i, P:-slopeOrderedIndices[i][P:-ns[i]]), vminus_index} , i = 0 .. P:-r)};
                    self:-Sigma := taus_plus union taus union taus_minus;
                end if;

            elif P:-classGroupRank = 1 then
                # In this case, the rank of the class group is one, hence there is only one complete fan in the lattice 
                # containing the columns of the P-Matrix as rays.
                self:-Sigma := {seq({seq(1 .. P:-numCols)} minus {i}, i = 1 .. P:-numCols)};
            elif admitsFano(P) then
                # In this case, we call the procedure again with the anticanonical class as the weight.
                return ComplexityOneVariety[ModuleCopy](self, proto, P, -getCanonicalDivisorClass(P));
            else
                error "This PMatrix is neither of Picard number one, nor is it a surface, nor does it admit a Fano variety."
                      "Therefore, you must provide a fan Sigma or a weight w as input.";
            end if;
        end if;
    end;

    export setMaximalXCones :: static := proc(self :: ComplexityOneVariety, maximalXCones :: set(set(integer))) 
        self:-maximalXCones := maximalXCones;
    end;

    export getMaximalXCones :: static := proc(self :: ComplexityOneVariety)
        if type(self:-maximalXCones, undefined) or 'forceCompute' in [_passed] then
            setMaximalXCones(self, getMaximalXConesFormat(self:-P:-format, self:-Sigma));
        end if;
        return self:-maximalXCones;
    end;

    export setMinimalAmbientFan :: static := proc(self :: ComplexityOneVariety, minimalAmbientFan :: FAN) 
        self:-minimalAmbientFan := minimalAmbientFan;
    end;

    export getMinimalAmbientFan :: static := proc(self :: ComplexityOneVariety)
        if type(self:-minimalAmbientFan, undefined) or 'forceCompute' in [_passed] then
            setMinimalAmbientFan(self, fan(op(map(c -> intSetConeToConvexCone(self:-P, c), getMaximalXCones(self)))));
        end if;
        return self:-minimalAmbientFan;
    end;

    export setIsQfactorialVal :: static := proc(self :: ComplexityOneVariety, isQfactorialVal :: boolean)
        self:-isQfactorialVal := isQfactorialVal
    end proc;

    export isQfactorial :: static := proc(self :: ComplexityOneVariety)
        if type(self:-isQfactorialVal, undefined) or 'forceCompute' in [_passed] then
            setIsQfactorialVal(self, issimplicial(getMinimalAmbientFan(self)));
        end if;
        return self:-isQfactorialVal;
    end proc;

    export setIsFactorialVal :: static := proc(self :: ComplexityOneVariety, isFactorialVal :: boolean)
        self:-isFactorialVal := isFactorialVal
    end proc;

    export isFactorial :: static := proc(self :: ComplexityOneVariety)
        if type(self:-isFactorialVal, undefined) or 'forceCompute' in [_passed] then
            setIsFactorialVal(self, isregular(getMinimalAmbientFan(self)));
        end if;
        return self:-isFactorialVal;
    end proc;

    export setLocalPicardIndices :: static := proc(self :: ComplexityOneVariety, localPicardIndices)
        self:-localPicardIndices := localPicardIndices;
    end proc;

    export getLocalPicardIndices :: static := proc(self :: ComplexityOneVariety)
        if type(self:-localPicardIndices, undefined) or 'forceCompute' in [_passed] then
            setLocalPicardIndices(self, map(c -> localPicardIndex(self:-P, c), convert(getMaximalXCones(self), list)));
        end if;
        return self:-localPicardIndices;
    end proc;

    export setPicardIndex :: static := proc(self :: ComplexityOneVariety, picardIndex)
        self:-picardIndex := picardIndex;
    end proc;

    export getPicardIndex :: static := proc(self :: ComplexityOneVariety)
        local P, Q, cones, clsGroup, localPicardGroups, picardGroup, factGrp, r, numTorsion, torsionColumns, lattices;
        
        if type(self:-picardIndex, undefined) or 'forceCompute' in [_passed] then
            # The picard index is the index of the picard group in the class group.
            # Note that it is *not* equal to the lcm of the local picard indices, since
            # the index of the intersection of two subgroups in an abelian group is not
            # equal to the lcm of the individual indices (the lcm is a lower bound though)

            r := getClassGroup(self:-P)[1];
            numTorsion := nops(getClassGroup(self:-P)) - 1;
            torsionColumns := <Matrix(r, numTorsion, fill = 0) ; Matrix(numTorsion, numTorsion, Vector(getClassGroup(self:-P)[2 .. -1]), shape = diagonal)>;

            cones := map(c -> {seq(1 .. self:-P:-numCols)} minus c, convert(getMaximalXCones(self), list));
            Q := getDegreeMatrix(self:-P);
            lattices := map(c -> <Matrix([Column(Q, convert(c, list))]) | torsionColumns>, cones);

            setPicardIndex(self, indexOfImage(integerIntersectionBasisList(lattices)));
        end if;
        return self:-picardIndex;
    end proc;

    export getLocalCartierIndices :: static := proc(self :: ComplexityOneVariety, D :: list(integer))
        map(c -> localCartierIndex(self:-P, c, D), convert(getMaximalXCones(self), list));
    end proc;

    export getCartierIndex :: static := proc(self :: ComplexityOneVariety, D :: list(integer))
        lcm(op(getLocalCartierIndices(self, D)));
    end proc;

    export setLocalGorensteinIndices :: static := proc(self :: ComplexityOneVariety, localGorensteinIndices :: list)
        self:-localGorensteinIndices := localGorensteinIndices;
    end proc;

    export getLocalGorensteinIndices :: static := proc(self :: ComplexityOneVariety)
        if type(self:-localGorensteinIndices, undefined) or 'forceCompute' in [_passed] then
            setLocalGorensteinIndices(self, map(c -> localGorensteinIndex(self:-P, c), convert(getMaximalXCones(self), list)));
        end if;
        return self:-localGorensteinIndices;
    end proc;

    export setGorensteinIndex :: static := proc(self :: ComplexityOneVariety, gorensteinIndex) 
        self:-gorensteinIndex := gorensteinIndex;
    end proc;

    export getGorensteinIndex :: static := proc(self :: ComplexityOneVariety)
        local cone;
        if type(self:-gorensteinIndex, undefined) or 'forceCompute' in [_passed] then
            setGorensteinIndex(self, lcm(op(getLocalGorensteinIndices(self))));
        end if;
        return self:-gorensteinIndex;
    end proc;

    export setIsGorensteinVal :: static := proc(self :: ComplexityOneVariety, isGorensteinVal :: boolean) 
        self:-isGorensteinVal := isGorensteinVal;
    end proc;

    export isGorenstein :: static := proc(self :: ComplexityOneVariety)
        local cone;
        if type(self:-isGorensteinVal, undefined) or 'forceCompute' in [_passed] then
            setIsGorensteinVal(self, evalb(getGorensteinIndex(self) = 1));
        end if;
        return self:-isGorensteinVal;
    end proc;

    export setIsQgorensteinVal :: static := proc(self :: ComplexityOneVariety, isQgorensteinVal :: boolean)
        self:-isQgorensteinVal := isQgorensteinVal;
    end proc;

    export isQgorenstein :: static := proc(self :: ComplexityOneVariety)
        if type(self:-isQgorensteinVal, undefined) or 'forceCompute' in [_passed] then
            setIsQgorensteinVal(self, evalb(getGorensteinIndex(self) < infinity));
        end if;
        return self:-isQgorensteinVal;
    end proc;

    export setLocalGorenstienQuotients :: static := proc(self :: ComplexityOneVariety, localGorensteinQuotients)
        self:-localGorensteinQuotients := localGorensteinQuotients;
    end proc;

    export getLocalGorensteinQuotients :: static := proc(self :: ComplexityOneVariety)
        if type(self:-localGorensteinQuotients, undefined) or 'forceCompute' in [_passed] then
            setLocalGorenstienQuotients(self, zip((x,y) -> x / y, getLocalPicardIndices(self), getLocalGorensteinIndices(self)));
        end if;
        return self:-localGorensteinQuotients;
    end proc;

    export setGorensteinQuotient :: static := proc(self :: ComplexityOneVariety, gorensteinQuotient)
        self:-gorensteinQuotient := gorensteinQuotient;
    end proc;

    export getGorensteinQuotient :: static := proc(self :: ComplexityOneVariety)
        local cone;
        if type(self:-gorensteinQuotient, undefined) or 'forceCompute' in [_passed] then
            setGorensteinQuotient(self, lcm(op(getLocalGorensteinQuotients(self))));
        end if;
        return self:-gorensteinQuotient;
    end proc;

    export setAmpleCone :: static := proc(self :: ComplexityOneVariety, ampleCone :: CONE) self:-ampleCone := ampleCone; end proc;

    export getAmpleCone :: static := proc(self :: ComplexityOneVariety)
        local cone;
        if type(self:-ampleCone, undefined) or 'forceCompute' in [_passed] then
            setAmpleCone(self, intersection(seq(poshull(Column(getDegreeMatrixFree(self:-P), [op({seq(1 .. ColumnDimension(getDegreeMatrixFree(self:-P)))} minus cone)])), cone in getMaximalXCones(self))));
        end if;
        return self:-ampleCone;
    end proc;

    export setIsFanoVal :: static := proc(self :: ComplexityOneVariety, isFanoVal :: boolean) self:-isFanoVal := isFanoVal; end proc;

    export isFano :: static := proc(self :: ComplexityOneVariety)
        if type(self:-isFanoVal, undefined) or 'forceCompute' in [_passed] then
            setIsFanoVal(self, containsrelint(getAmpleCone(self), -getCanonicalDivisorClass(self:-P)));
        end if;
        return self:-isFanoVal
    end proc;

    export setIntersectionMatrix :: static := proc(self :: ComplexityOneVariety, intersectionMatrix :: Array) self:-intersectionMatrix := intersectionMatrix; end proc;

    export getIntersectionMatrix :: static := proc(X :: ComplexityOneVariety)   
        local P, res, i, mcal, newM, j, j_, j1, j2, k1, k2, k, i1, i2, kplus, kminus;

        if type(X:-intersectionMatrix, undefined) or 'forceCompute' in [_passed] then
            P := X:-P;
            if X:-P:-s > 1 then
                error "Intersection numbers are only defined for K*-surfaces, i.e. s = 1.";
            end if;

            # We use the formulas of Chapter 6 in the paper: "Log del pezzo surfaces with torus action" by Hättig, Hausen, Hummel
            # Note however that we do not require P to be slope-ordered, hence we have to generalize the formulas
            # a bit by translating the indices through the `P:-slopeOrderedIndices` array.

            res := Array(1 .. P:-numCols, 1 .. P:-numCols, fill = 0);
            P:-slopeOrderedIndices := Array(0..P:-r, [seq(sort([seq(1 .. P:-ns[i])], 
                        (j1, j2) -> P:-slopes[i][j1] > P:-slopes[i][j2]), 
                        i = 0 .. P:-r)]);
            
            # First, compute calligraphic m. This is encoded as an array of arrays, since we want to have
            # the indexing go from 0 to n_i.
            mcal := Array(0 .. P:-r);
            for i from 0 to P:-r do
                mcal[i] := Array(0 .. P:-ns[i], fill = 0);
                if P:-case = "EE" or P:-case = "EP" then
                    # There is an elliptic fixed point x^+
                    mcal[i][0] := - 1 / P:-mplus;
                end if;
                if P:-case = "EE" or P:-case = "PE" then
                    # There is an elliptic fixed point x^-
                    # (Different from printed formula, since our mminus is negated)
                    mcal[i][P:-ns[i]] := - 1 / P:-mminus;
                end if;

                for j from 1 to P:-ns[i] - 1 do
                    mcal[i][j] := 1 / (P:-slopes[i][P:-slopeOrderedIndices[i][j]] - P:-slopes[i][P:-slopeOrderedIndices[i][j+1]]);
                end do;
            end do;

            # Compute intersection numbers of two adjacent rays in the leaves.
            # This is independent of the case of P
            for i from 0 to P:-r do
                for j_ from 1 to P:-ns[i] - 1 do
                    j1 := P:-slopeOrderedIndices[i][j_];
                    j2 := P:-slopeOrderedIndices[i][j_+1];
                    k1 := doubleToSingleIndex(P:-format, i, j1);
                    k2 := doubleToSingleIndex(P:-format, i, j2);
                    res[k1,k2] := 1 / (P:-lss[i][j1] * P:-lss[i][j2]) * mcal[i][j_];
                    res[k2,k1] := res[k1,k2];
                end do;
            end do:

            # Compute the self intersection numbers of rays in the leaves
            # This is independent of the case of P
            for i from 0 to P:-r do
                for j_ from 1 to P:-ns[i] do
                    j := P:-slopeOrderedIndices[i][j_];
                    k := doubleToSingleIndex(P:-format, i, j);
                    res[k,k] := - 1 / P:-lss[i][j]^2 * (mcal[i][j_ - 1] + mcal[i][j_]);
                end do;
            end do:

            if P:-case = "EE" or P:-case = "EP" then
                ##########################################  
                ## There is an elliptic fixed point x^+ ##
                ##########################################

                # Intersection numbers of the highest rays in each block with each other.
                for i1 from 0 to P:-r do
                    for i2 in {seq(0 .. P:-r)} minus {i1} do
                        j1 := P:-slopeOrderedIndices[i1][1];
                        j2 := P:-slopeOrderedIndices[i2][1];
                        k1 := doubleToSingleIndex(P:-format, i1, j1);
                        k2 := doubleToSingleIndex(P:-format, i2, j2);
                        if P:-ns[i1] = 1 and P:-ns[i2] = 1 then
                            res[k1,k2] := - 1 / (P:-lss[i1][j1] * P:-lss[i2][j2]) * (mcal[i1][0] + mcal[i1][1]);
                        else 
                            res[k1,k2] := - 1 / (P:-lss[i1][j1] * P:-lss[i2][j2]) * mcal[i1][0];
                        end if;
                    end do;
                end do;
            else
                ################################################
                ## There is a parabolic fixed point curce D^+ ##
                ################################################

                # kplus is the index of the divisor D^+.
                if P:-d[1, P:-n + 1] = 1 then
                    kplus := P:-n + 1;
                elif P:-d[1, P:-n + 1] = -1 then
                    kplus := P:-n + 2;
                end if;

                # Intersection numbers of the highest rays of each block with D^+
                for i from 0 to P:-r do
                    j := P:-slopeOrderedIndices[i][1];
                    k := doubleToSingleIndex(P:-format, i, j);
                    res[k,kplus] := 1 / P:-lss[i][j];
                    res[kplus,k] := res[k,kplus];
                end do;

                # Self intersection of D^+
                res[kplus, kplus] := - P:-mplus;

            end if;

            if P:-case = "EE" or P:-case = "PE" then
                ##########################################  
                ## There is an elliptic fixed point x^- ##
                ##########################################

                # Intersection numbers of the lowest rays in each block with each other.
                for i1 from 0 to P:-r do
                    for i2 in {seq(0 .. P:-r)} minus {i1} do
                        j1 := P:-slopeOrderedIndices[i1][P:-ns[i1]];
                        j2 := P:-slopeOrderedIndices[i2][P:-ns[i2]];
                        k1 := doubleToSingleIndex(P:-format, i1, j1);
                        k2 := doubleToSingleIndex(P:-format, i2, j2);
                        if P:-ns[i1] = 1 and P:-ns[i2] = 1 then
                            res[k1,k2] := - 1 / (P:-lss[i1][j1] * P:-lss[i2][j2]) * (mcal[i1][0] + mcal[i1][1]);
                        else 
                            res[k1,k2] := - 1 / (P:-lss[i1][j1] * P:-lss[i2][j2]) * mcal[i1][P:-ns[i1]];
                        end if;
                    end do;
                end do;
            else 
                ################################################
                ## There is a parabolic fixed point curce D^- ##
                ################################################

                # kminus is the index of the divisor D^-.
                if P:-d[1, P:-n + 1] = -1 then
                    kminus := P:-n + 1;
                elif P:-d[1, P:-n + 1] = 1 then
                    kminus := P:-n + 2;
                end if;

                # Intersection numbers of the highest rays of each block with D^-
                for i from 0 to P:-r do
                    j := P:-slopeOrderedIndices[i][P:-ns[i]];
                    k := doubleToSingleIndex(P:-format, i, j);
                    res[k,kminus] := 1 / P:-lss[i][j];
                    res[kminus,k] := res[k,kminus];
                end do;

                # Self intersection of D^-
                # (Different from printed formula since out mminus is negated)
                res[kminus, kminus] := - P:-mminus;

            end if;

            setIntersectionMatrix(X, res);
        
        end if;

        return X:-intersectionMatrix;

    end proc;

    (*
    Compute the intersection number of any two divisors, given as linear combinations
    of the primitive invariant divisors D_X^{ij} and D_X^{\pm}.
    The input data are lists of integers encoding the coefficients this linear combination.
    *)
    export intersectionNumber :: static := proc(X :: ComplexityOneVariety, D1 :: list(integer), D2 :: list(integer))
        local k1, k2;
        
        if nops(D1) <> X:-P:-numCols or nops(D2) <> X:-P:-numCols then
            error "The list of integers encoding the divisor must have length n + m = %1", X:-P:-n + X:-P:-m;
        end if;

        return add([seq(seq(D1[k1] * D2[k2] * getIntersectionMatrix(X)[k1,k2], k2 = 1 .. nops(D2)), k1 = 1 .. nops(D1))]);

    end proc;

    export setAnticanonicalSelfIntersection :: static := proc(self :: ComplexityOneVariety, anticanonicalSelfIntersection) self:-anticanonicalSelfIntersection := anticanonicalSelfIntersection; end proc;

    export getAnticanonicalSelfIntersection :: static := proc(X :: ComplexityOneVariety)
        if type(X:-anticanonicalSelfIntersection, undefined) or 'forceCompute' in [_passed] then
            setAnticanonicalSelfIntersection(X, intersectionNumber(X, -getCanonicalDivisorCoefficients(X:-P), -getCanonicalDivisorCoefficients(X:-P)));
        end if;
        return X:-anticanonicalSelfIntersection;
    end proc;


    (****************************
    *** ADMISSIBLE OPERATIONS ***
    *****************************)

    export removeSingleErasableBlock :: static := proc(X :: ComplexityOneVariety, i0 :: integer)
        local newP, oldToNewIndex, newSigma;
        newP := PMatrix[removeSingleErasableBlock](X:-P, i0);
        oldToNewIndex := k -> if k < add(X:-P:-ns[0 .. i0 - 1]) + 1 then k else k - 1 end if;
        newSigma := map(cones -> map(oldToNewIndex, cones), X:-Sigma);
        return ComplexityOneVariety(newP, newSigma);
    end proc;

    export removeErasableBlocks :: static := proc(X :: ComplexityOneVariety)
        local P, redundantIndices;
        P := X:-P;
        redundantIndices := select(i -> P:-lss[i] = [1], [seq(0 .. P:-r)]);
        # If there is a redundant block and we still have more than two blocks, remove it.
        if nops(redundantIndices) > 0 and P:-r > 1 then
            # Recursive call
            return removeErasableBlocks(removeSingleErasableBlock(X, redundantIndices[1]));
        else
            return X;
        end if;
    end proc;

    export applyAdmissibleOperation :: static := proc(X :: ComplexityOneVariety, a :: AdmissibleOperation)
        local newP, newSigma, newA;
        newP := PMatrix[applyAdmissibleOperation](X:-P, a);
        newSigma := map(cones -> map(k -> a:-bundledPermutation[k], cones), X:-Sigma);
        if type(X:-A, undefined) then
            ComplexityOneVariety(newP, newSigma);
        else
            newA := a:-U . X:-A . a:-D . a:-sigmaPermutationMatrix^(-1);
            ComplexityOneVariety(newP, newSigma, newA);
        end if;

    end proc;

    export sortColumnsByLss :: static := proc(X :: ComplexityOneVariety)
       applyAdmissibleOperation(X, PMatrix[sortColumnsByLss](X:-P, 'operation'));
    end proc;

    export standardizeCoefficientMatrixOperation :: static := proc(X :: ComplexityOneVariety)
        local U1, newA, ds, U2, i;
        U1 := Matrix([Column(X:-A,[1,2])])^(-1);
        newA := U1 . X:-A;
        ds := [1, newA[2,3] / newA[1,3], seq(-1 / newA[1,i], i = 3 .. X:-P:-r + 1)];
        U2 := <1,0;0,newA[1,3]/newA[2,3]>;
        AdmissibleOperation[OnA](X:-P:-format, U2 . U1, ds);
    end proc;

    export standardizeCoefficientMatrix :: static := proc(X :: ComplexityOneVariety)
        applyAdmissibleOperation(X, standardizeCoefficientMatrixOperation(X));
    end proc;
    
    export areCoefficientMatricesEquivalent :: static := proc(X1 :: ComplexityOneVariety, X2 :: ComplexityOneVariety)
        local a1, a2, resultOp, out;

        # Output handling
        if _npassed > 2 then
            if not (_passed[3] in {'boolean', 'operation'}) then
                error "third parameter must be either 'boolean' or 'operation'";
            end if;
            out := _passed[3]
        else
            out := 'boolean';
        end if;

        a1 := standardizeCoefficientMatrixOperation(X1);
        a2 := standardizeCoefficientMatrixOperation(X2);
        resultOp := NULL;
        if Equal(applyAdmissibleOperation(X1, a1):-A, applyAdmissibleOperation(X2, a2):-A) then
            resultOp := compose(a1, inverse(a2));
        end if;
       
        if out = 'operation' then
            return resultOp;
        else
            return type(resultOp, AdmissibleOperation);
        end if;
    end proc;

    (*
    Checks whether two complexity-one varieties are isomorphic to each other
    *)
    export areIsomorphic :: static := proc(X1_ :: ComplexityOneVariety, X2_ :: ComplexityOneVariety)
        local X1, X2, a, newX1, coefficientOp, out, resultOps;

        # Output handling
        if _npassed > 2 then
            if not (_passed[3] in {'boolean', 'operation'}) then
                error "third parameter must be either 'boolean' or 'operation'";
            end if;
            out := _passed[3]
        else
            out := 'boolean';
        end if;

        X1 := removeErasableBlocks(X1_);
        X2 := removeErasableBlocks(X2_);

        resultOps := [];

        for a in PMatrix[areEquivalent](X1:-P, X2:-P, 'operation') do
            newX1 := applyAdmissibleOperation(X1, a);
            if getMaximalXCones(newX1) = getMaximalXCones(X2) then
                # If there are no coefficient matrices, we are done.
                if type(X1:-A, undefined) or type(X2:-A, undefined) then
                    if out = 'boolean' then
                        return true;
                    end if;
                    resultOps := [op(resultOps), a];
                else
                    # Else we check if the coefficient matrices can be turned into each other by an
                    # admissible operation on the coefficient matrix.
                    coefficientOp := areCoefficientMatricesEquivalent(newX1, X2, 'operation');
                    if type(coefficientOp, AdmissibleOperation) then
                        if out = 'boolean' then
                            return true;
                        end if;
                        resultOps := [op(resultOps), compose(a, coefficientOp)];
                    end if;
                end if;
            end if;
        end do;

        if out = 'operation' then
            return resultOps;
        else
            return evalb(resultOps <> []);
        end if;

    end proc;

    export normalForm :: static := proc(X :: ComplexityOneVariety)
        applyAdmissibleOperation(X, PMatrix[normalForm](X:-P, 'operation'));
    end proc;

    (*
    Construct the tropical resolution of a ComplexityOneVariety.
    *)
    export tropicalResolution :: static := proc(X :: ComplexityOneVariety)
        local P, tropicalSheetsGenerators, tropicalCones, newCones, Pcols, newPMatrixColumns, newP, newSigma, k, i;

        P := X:-P;

        # The tropical cones are the cones rho_i x QQ^s, where rho_i is generated by the canonical basis vector e_i
        # and e_0 = -(e_1 + ... + e_r). 
        tropicalSheetsGenerators := [seq([0 $ (P:-r + k - 1), 1, 0 $ (P:-s - k)], k = 1 .. P:-s), seq([0 $ (P:-r + k - 1), -1, 0 $ (P:-s - k)], k = 1 .. P:-s)];
        tropicalCones := [seq(poshull(convert(<canonicalBasisVector(P:-r, i) ; Vector(P:-s, fill = 0)>, list), op(tropicalSheetsGenerators)), i = 0 .. P:-r)];

        # To get the new cones of our fan, we intersect each X-Cone with each of the tropical cones.
        # This might yield new rays that we append to our P-Matrix later
        newCones := map(c1 -> op(map(c2 -> intersection(c1,c2), tropicalCones)), map(c -> intSetConeToConvexCone(X:-P, c), getMaximalXCones(X)));
        
        # We construct newP, which arises from P by appending all the columns that arised during the intersection with the tropical cones.
        Pcols := map(v -> :-convert(v, list), [Column(P:-mat, [seq(1 .. P:-n + P:-m)])]);
        newPMatrixColumns := select(ray -> not (ray in Pcols), map(c -> op(rays(c)), newCones));
        newP := PMatrix(P:-s, <P:-mat | Transpose(Matrix([op(newPMatrixColumns)]))>);

        # Convert back from CONEs to lists of integers to get out new fan
        newSigma := map(c -> convexConeToIntSetCone(newP, c), newCones);

        return ComplexityOneVariety(newP, newSigma);

    end proc;

    export canonicalResolution :: static := proc(X0 :: ComplexityOneVariety)
        local X, P, newFan, Pcols, allRays, newRaysByBlocks, leaf, newPColumns, newP, newSigma, newX, i, exceptDivIndices;

        X := tropicalResolution(X0);
        P := X:-P;
        newFan := regularsubdiv(getMinimalAmbientFan(X));

        Pcols := map(v -> :-convert(v, list), [Column(P:-mat, [seq(1 .. P:-numCols)])]);
        allRays := {op(map(c -> op(rays(c)), maximal(newFan)))};

        newRaysByBlocks := Array(0..P:-r);
        for i from 0 to P:-r do
            leaf := poshull(convert(canonicalBasisVector(P:-r, i), list));
            newRaysByBlocks[i] := select(ray -> not (ray in Pcols) and containsrelint(leaf, ray[1 .. P:-r]), allRays);
        end do;

        newPColumns := [];
        exceptDivIndices := [];
        for i from 0 to P:-r do
            # First add the old P-columns
            newPColumns := [op(newPColumns), op(map(v -> :-convert(v, list), [Column(P:-mat, [seq(add(P:-ns[0 .. i-1]) + 1 .. add(P:-ns[0 .. i]))])]))];
            exceptDivIndices := [op(exceptDivIndices), seq(nops(newPColumns) + 1  .. nops(newPColumns) + nops(newRaysByBlocks[i]))];
            # now append the new ones at the end of the block
            newPColumns := [op(newPColumns), op(newRaysByBlocks[i])];
        end do;
        newPColumns := [op(newPColumns), op(map(v -> :-convert(v, list), [Column(P:-mat, [seq(P:-n + 1 .. P:-numCols)])]))];

        newP := PMatrix(P:-s, Transpose(Matrix(newPColumns)));
        newSigma := {op(map(c -> convexConeToIntSetCone(newP, c), maximal(newFan)))};
        newX := ComplexityOneVariety(newP, newSigma);
        newX:-exceptionalDivisorsIndices := exceptDivIndices;
        return newX;

    end proc;

    export minimalResolution :: static := proc(X0 :: ComplexityOneVariety)
        local X, P, contractibleIndices;
        X := canonicalResolution(X0);
        contractibleIndices := select(i -> getIntersectionMatrix(X)[i,i] = -1, X:-exceptionalDivisorsIndices);
        while contractibleIndices <> [] do
            X := ComplexityOneVariety(PMatrix(1, DeleteColumn(X:-P:-mat, contractibleIndices)));
            contractibleIndices := select(i -> getIntersectionMatrix(X)[i,i] = -1, X:-exceptionalDivisorsIndices);
        end do;
        return X;
    end proc;

    export ModulePrint :: static := proc(self :: ComplexityOneVariety)
        nprintf(cat("ComplexityOneVariety(dim = ", self:-P:-s + 1,
          ", lss = ", convert(self:-P:-lss, list),
          ", Sigma = ", self:-Sigma));
    end;

    export ComplexityOneVarietyInfo :: static := proc(self :: ComplexityOneVariety)
        local P, i, relations, maximalXCones, Q, classGroup, classGroupRank, anticanonicalClass, effectiveConeRays, movingConeRays, ampleConeRays, isFano, gorensteinIndex, intersectionMatrix, anticanonicalSelfIntersection;
        print(P = self:-P:-mat);
        print([seq(cat(n,i), i = 0 .. self:-P:-r - 1), m] = [seq(self:-P:-ns[i], i = 1 .. self:-P:-r), self:-P:-m]);
        print(relations = self:-relations);
        print(maximalXCones = getMaximalXCones(self));
        print(Q = getDegreeMatrix(self:-P));
        print(classGroup = getClassGroup(self:-P));
        print(classGroupRank = self:-P:-classGroupRank);
        print(canonicalClass = getCanonicalDivisorClass(self:-P));
        print(effectiveConeRays = rays(getEffectiveCone(self:-P)));
        print(movingConeRays = rays(getMovingCone(self:-P)));
        print(ampleConeRays = rays(getAmpleCone(self)));
        print(isFano = ComplexityOneVariety[isFano](self));
        print(gorensteinIndex = getGorensteinIndex(self));
        if self:-P:-s = 1 then
            print(intersectionMatrix = getIntersectionMatrix(self));
            print(anticanonicalSelfIntersection = getAnticanonicalSelfIntersection(self));
        end if;
    end;

end module: