with(PDEtools): with(DEtools):

## Hirota Bilinear Method
## Bilinear Derivative / Hirota Operator
BD := proc(FF, DD)
    local f, g, x, m, opt;
    if nargs = 1 then
        return `*`(FF[]);
    fi;
    f, g := FF[];
    x, m := DD[];
    opt := args[3..-1];
    if m = 0 then
        return procname(FF, opt);
    fi;
    procname([diff(f, x), g], [x, m-1], opt) - procname([f, diff(g, x)], [x, m-1], opt);
end:

`print/BD` := proc(FF, DD)
    local f, g, x, m, i;
    f, g := FF[];
    f := cat(f, ` `, g);
    g := product(D[args[i][1]]^args[i][2], i = 2..nargs);
    if g <> 1 then
        f := ``(g)*``(f);
    fi;
    f;
end:

## collect(expr, f); first!
getFnumer := proc(df, f, pow::posint := 1)
    local i, g, fdenom;
    if type(df, `+`) then
        g := [op(df)];
        fdenom := map(denom, g);
        for i to nops(fdenom) while fdenom[i] <> f^pow do
        od;
        if i > nops(fdenom) then
            lprint(fdenom);
            error "no term(s) or numer=0 when denom=%1", op(0, f)^pow;
        fi;
        g := numer(g[i]);
        if not type(expand(g), `+`) then
            lprint(g);
            error "Expected more than 1 term about Hirota D-operator";
        fi;
        return g;
    fi;
    lprint(df);
    error "expected 1st argument be type `+`.";
end:

getvarpow := proc(df::function)
    local i, f, var, dif, pow;
    if op(0, df) <> diff then
        lprint(df);
        error "expected diff function";
    fi;
    f := convert(df, D);
    var := [op(f)];
    dif := [op(op([0, 0], f))];
    pow := [0$nops(var)];
    f := op(op(0, f))(var[]);
    for i to nops(var) do
        dif := selectremove(member, dif, {i});
        pow[i] := nops(dif[1]);
        dif := dif[2];
    od;
    pow := zip((x, y) -> [x, y], var, pow);
    pow := remove(has, pow, {0});
    [[f, f], pow[]];
end:

#convert to Hirota Bilinear Form
HBF := proc(df)
    local i, c, f;
    if type(df, `+`) then
        f := [op(df)];
        return map(procname, f);
    fi;
    if type(df, `*`) then
        f := [op(df)];
        f := selectremove(hasfun, f, diff);
        c := f[2];
        f := f[1];
        if nops(f) <> 1 then
            lprint(df);
            error "need only one diff function factor.";
        fi;
        f := f[];
        c := `*`(c[]);
        f := getvarpow(f);
        f := [c, f];
        return f;
    fi;
    if op(0, df) = diff then
        f := getvarpow(df);
        f := [1, f];
        return f;
    fi;
    lprint(df);
    error "unexpected type.";
end:

printHBF := proc(PL::list)
    local j, DD, f, C, tmp, gcdC, i;
    C := map2(op, 1, PL);
    gcdC := 1;
    if nops(C) > 1 then
        tmp := [seq(cat(_Z, i), i = 1..nops(C))];
        gcdC := tmp *~ C;
        gcdC := `+`(gcdC[]);
        gcdC := factor(gcdC);
        tmp := selectremove(has, gcdC, tmp);
        gcdC := tmp[2];
        if gcdC = 0 then
            gcdC := 1;
        fi;
        gcdC := gcdC*content(tmp[1]);
    fi;
    if gcdC <> 1 then
        C := C /~ gcdC;
    fi;
    DD := map2(op, 2, PL);
    f := op(0, DD[1][1][1]);
    DD := map(z -> product(D[z[i][1]]^z[i][2], i = 2..nops(z)), DD);
    DD := zip(`*`, C, DD);
    DD := `+`(DD[]);
    gcdC * ``(DD) * cat(f, ` `, f);
end:

## print Hirota Bilinear Transform
printHBT := proc(uf, u, f, i, j, PL, alpha := 1)
    local DD, g, C, tmp, pl;
    pl := printHBF(PL);
    if j > 0 then
        print(u = 2*alpha*'diff'(ln(f), x$j));
    else
        print(u = 2*alpha*ln(f));
    fi;
    if i > 0 then
        print('diff'(pl/f^2, x$i) = 0);
    else
       print(pl/f^2 = 0);
    fi;
    NULL;
 end:

guessdifforder := proc(PL::list, x::name)
    local L, minorder, maxorder, tmp;
    L := map2(op, 2, PL);
    L := map(z -> z[2..-1], L);
    tmp := map(z -> map2(op, 2, z), L);
    tmp := map(z -> `+`(z[]), tmp);
    tmp := selectremove(type, tmp, even);
    minorder := 0;
    if nops(tmp[1]) < nops(tmp[2]) then
        minorder := 1;
    fi;
    tmp := map(z -> select(has, z, {x}), L);
    tmp := map(z -> map2(op, 2, z), tmp);
    if has(tmp, {[]}) then
        maxorder := 0;
    else
        tmp := map(op, tmp);
        maxorder := min(tmp[]);
    fi;
    if type(maxorder - minorder, odd) then
        maxorder := maxorder - 1;
    fi;
    [minorder, maxorder];
end:

guessalpha := proc(Res, uf, u, f, i, j, PL, alpha)
    local tmp, res, pl, flag, k;
    flag := 1;
    tmp := [op(Res)];
    tmp := map(numer, tmp);
    tmp := gcd(tmp[1], tmp[-1]);
    if type(tmp, `*`) then
        tmp := remove(has, tmp, f);
    fi;
    if tmp <> 0 and has(tmp, {alpha}) then
        tmp := solve(tmp/alpha^difforder(uf), {alpha});
        if tmp <> NULL and has(tmp, {alpha}) then
            lprint(tmp);
            for k to nops([tmp]) while flag = 1 do
                res := collect(expand(subs(tmp[k], Res)), f, factor);
                if res = 0 then
                    pl := subs(tmp[k], PL);
                    printHBT(uf, u, f, i, j, pl, rhs(tmp[k]));
                    flag := 0;
                fi;
            od;
        fi;
    fi;
    PL;
end:

Bilinear := proc(uf, u, f, x, alpha)
    local su, h, i, j, g1, CB, PL, gdo, DD, Res;
    if hasfun(uf, int) then
        error "Do not support integral function yet. Please substitute int function.";
    fi;
    for j from 0 to 2 do
        Res := 1;
        su := u = 2*alpha*diff(ln(f), [x$j]);
        h := collect(expand(dsubs(su, uf)), f, factor);
        if hasfun(h, ln) then
            next;
        fi;
        g1 := getFnumer(h, f)/2;
        g1 := expand(g1);
        CB := HBF(g1);
        gdo := guessdifforder(CB, x);
        for i from gdo[1] by 2 to gdo[2] do
            if i = 0 then
                PL := CB;
            else
                PL := HBF(int(g1, x$i));
            fi;
            DD := add(PL[i][1]*BD(PL[i][2][]), i = 1..nops(PL));
            Res := collect(expand(diff(DD/f^2, [x$i]) - h), f, factor);
            if Res = 0 then
                printHBT(uf, u, f, i, j, PL, alpha);
                break;
            elif type(alpha, name) and has(DD, alpha) then
                Res := guessalpha(Res, uf, u, f, i, j, PL, alpha);
            fi;
        od;
        if Res = 0 then
            break;
        fi;
    od;
    PL;
end:

with(PDEtools):

# Definieer de Boussinesq vergelijking
boussinesq := diff(u(x,t),t,t) - diff(u(x,t),x,x) - 3*diff(u(x,t),x,x)^2 - diff(u(x,t),x$4);

# Pas de bilineaire transformatie toe
Bilinear(boussinesq, u(x,t), f(x,t), x, alpha);

{alpha = 1}

# Definieer de Kadomtsev–Petviashvili (KP) vergelijking in één regel
kp := diff(diff(u(x,y,t),x$3) + 6*u(x,y,t)*diff(u(x,y,t),x) + diff(u(x,y,t),t), x)
       + 3*delta^2*diff(u(x,y,t), y$2):
Bilinear(kp, u(x,y,t), f(x,y,t), x, alpha);
 

{alpha = 1}