kaleidicassociates / lubeck

High level linear algebra library for Dlang
http://lubeck.libmir.org/
Boost Software License 1.0
66 stars 14 forks source link

add pinv to lubeck2 and make svd nothrow #55

Open aferust opened 3 months ago

aferust commented 3 months ago

pinv is added to lubeck2. for potential code reviewers, please confirm this:

instead of just enforce!("svd: " ~ msg)(!info);

I had to use the below version to make the function nothrow. I am using this method in DCV a lot. It might be a problem for release builds due to assert, but at least it is nothrow, and works with debug builds.

try enforce!("svd: " ~ msg)(!info); catch(Exception e) assert(false, e.msg);

aferust commented 3 months ago

@9il please review this.

9il commented 3 months ago

What is wrong with mldivide?

9il commented 3 months ago

What is wrong with mldivide?

Ah, I see, mldivide will require identity matrix allocation. OK then

aferust commented 3 months ago

Hello, @9il, I am adapting an image stitching method (ransac homography estimation) DCV. I want to keep nogc nothrow with DCV. So, if we cannot have a function pinv nogc nothrow, I cannot do that. I have already tested my modification and I can do some elegant panaroma stitching. But yes I will add a unittest. If we don't have a nothrow solution, I will use my fork as a dependency?

9il commented 3 months ago

You didn't understand me. I suggest you to add an additional low-level API that will be ~throw~ nothrow and will incorporate the info, rework pinv based on this API, call this API inside the old API. It is a simple nothrow.

9il commented 3 months ago

damn, autocorrection replaces nothrow with throw

aferust commented 3 months ago

Ok, now I am working on the RANSAC thing. When I find some time, I will review the "info" and try to be more familiar with it and add an additoinal low-level API as you suggest. Until it, let's keep this PR open.

aferust commented 3 months ago

just to see what I was trying to do: https://github.com/aferust/dcv/tree/master/examples/imagestitchinghomography

aferust commented 3 months ago

You didn't understand me. I suggest you to add an additional low-level API that will be ~throw~ nothrow and will incorporate the info, rework pinv based on this API, call this API inside the old API. It is a simple nothrow.

do you mean creating a nothrow duplicate of svd and call it inside the pinv?

9il commented 3 months ago

You didn't understand me. I suggest you to add an additional low-level API that will be ~throw~ nothrow and will incorporate the info, rework pinv based on this API, call this API inside the old API. It is a simple nothrow.

do you mean creating a nothrow duplicate of svd and call it inside the pinv?

Yes. Please avoid code duplication. The old API should reuse the new one nothrow API.

aferust commented 3 months ago

I understood. I will make the required changes when i am available.

aferust commented 3 months ago

Before starting to do anything, I have just used a test code to determine an assertion workflow, but without touching anything, actually running on run.dlang.org the second assert fails. is it normal:

/+dub.sdl:
dependency "lubeck" version="~>1.5.4"
+/
import std.stdio;
import mir.ndslice, mir.rc;

void main()
{
    import kaleidic.lubeck;

    import std.math : isClose;
    auto A = iota(15).as!double.sliced(3,5) + 1;

    auto A_pinv = pinv(A);

    auto A1 = mtimes(A, mtimes(A_pinv, A));
    auto A2 = mtimes(A_pinv, mtimes(A, A_pinv));

    auto boolMat1 = zip(A,      A1).map!((a,b) => a.isClose(b));
    auto boolMat2 = zip(A_pinv, A2).map!((a,b) => a.isClose(b));

    // for potential asserts
    boolMat1.all!(aa => aa == true).writeln; // true
    boolMat2.all!(aa => aa == true).writeln; // false ? only one value is different

    A_pinv.writeln;    
    A2.writeln;
}

outputs: true false [[-0.246667, -0.0666667, 0.113333], [-0.133333, -0.0333333, 0.0666667], [-0.02, 3.77909e-18, 0.02], [0.0933333, 0.0333333, -0.0266667], [0.206667, 0.0666667, -0.0733333]] [[-0.246667, -0.0666667, 0.113333], [-0.133333, -0.0333333, 0.0666667], [-0.02, -2.54426e-18, 0.02], [0.0933333, 0.0333333, -0.0266667], [0.206667, 0.0666667, -0.0733333]]

jmh530 commented 3 months ago

Try this:

/+dub.sdl:
dependency "lubeck" version="~>1.5.4"
+/
import kaleidic.lubeck;

void main()
{
    import mir.algorithm.iteration: equal;
    import mir.math.common: approxEqual;
    import mir.ndslice.slice: sliced;
    import mir.ndslice.topology: iota, as;

    auto A = iota([15], 1).as!double.sliced(3,5);

    auto A_pinv = pinv(A);

    auto A1 = mtimes(A, mtimes(A_pinv, A));
    auto A2 = mtimes(A_pinv, mtimes(A, A_pinv));

    assert(A.equal!approxEqual(A1));
    assert(A_pinv.equal!approxEqual(A2));
}
aferust commented 3 months ago

@9il , @jmh530 , please review my last commit.

aferust commented 3 months ago

I understood, but I am outside now. I do it when İ get home

11 May 2024 Cmt 14:00 tarihinde Ilia Ki @.***> şunu yazdı:

@.**** commented on this pull request.

In source/kaleidic/lubeck2.d https://github.com/kaleidicassociates/lubeck/pull/55#discussion_r1597425426 :

@@ -1228,6 +1228,28 @@ Returns: error code from CBlas ) if (algorithm == "gesvd" || algorithm == "gesdd") { +

  • size_t info;
  • auto svdresult = svdImpl!(T, algorithm, kind)(a, info, slim);
  • enum msg = (algorithm == "gesvd" ? "svd: DBDSDC did not converge, updating process failed" : "svd: DBDSQR did not converge");
  • enforce!("svd: " ~ msg)(!info);
  • return svdresult; //transposed +}
  • @.*** pure @nogc nothrow SvdResult!T svdImpl

I think just svd name overload is OK.

— Reply to this email directly, view it on GitHub https://github.com/kaleidicassociates/lubeck/pull/55#pullrequestreview-2051181566, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACLDF6IGJM3NPUSJXCH3XK3ZBX25XAVCNFSM6AAAAABHQJZQR2VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDANJRGE4DCNJWGY . You are receiving this because you authored the thread.Message ID: @.***>

aferust commented 3 months ago

@9il I have added some pinv overloads. They cover possible situations. And there are no signature collisions.

9il commented 3 months ago

Folks, it seems I no longer have write access to Lubeck. End of story.

aferust commented 3 months ago

Folks, it seems I no longer have write access to Lubeck. End of story.

I am sorry to hear it. I hope things positively change over time.

jmh530 commented 3 months ago

@9il It's boost licensed. Would you want to fork it and do new work there?

9il commented 3 months ago

I prefer not to maintain the repository. I still have the code.dlang.org name lubeck.

First, we need to ask @Laeeth, if and where he wishes me to transfer the ownership of the lubeck package name on code.dlang.org. We will need to wait a while to get a response.

The safest thing for now is for @jmh530 or @aferust to create a fork with another package name and all credentials and links to this repository.

jmh530 commented 3 months ago

@9il I have a fork already. It doesn't let you fork it twice.