Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add pinv to lubeck2 and make svd nothrow #55

Open
wants to merge 7 commits into
base: master
Choose a base branch
from

Conversation

aferust
Copy link

@aferust aferust commented May 10, 2024

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
Copy link
Author

aferust commented May 10, 2024

@9il please review this.

enforce!("svd: " ~ msg)(!info);

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

@9il 9il May 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

SVD isn't supposed to be nothrow with the current API. You may want to add an additoinal low-level API that will incorporate the info into the result.

{
auto matrixScope = matrix.lightScope.lightConst;
return pinv(matrixScope, tolerance);
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please add a test

@9il
Copy link
Contributor

9il commented May 10, 2024

What is wrong with mldivide?

@9il
Copy link
Contributor

9il commented May 10, 2024

What is wrong with mldivide?

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

@aferust
Copy link
Author

aferust commented May 10, 2024

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
Copy link
Contributor

9il commented May 10, 2024

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
Copy link
Contributor

9il commented May 10, 2024

damn, autocorrection replaces nothrow with throw

@aferust
Copy link
Author

aferust commented May 10, 2024

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
Copy link
Author

aferust commented May 10, 2024

@aferust
Copy link
Author

aferust commented May 10, 2024

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?

else
auto si = st[0];
auto s = svd.sigma[0 .. $ - si];
s.each!"a = 1 / a";
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't tried this code or anything, but I would definitely want to make sure this line shouldn't be something like s.each!"a = 1.0 / a"

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

1 / is fine. It can be 1f / but 1.0 / . The last one will trigger float to double conversion.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree 1f is better. I was a little concerned because the function signature doesn't constrain T to be floating point (though I imagine some of the other functions called make that restriction, it isn't immediately obvious from looking at this one if that is the case).

@9il
Copy link
Contributor

9il commented May 10, 2024

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
Copy link
Author

aferust commented May 10, 2024

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

@aferust
Copy link
Author

aferust commented May 10, 2024

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
Copy link
Contributor

jmh530 commented May 10, 2024

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
Copy link
Author

aferust commented May 10, 2024

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

)(
auto ref const Slice!(const(T)*, 2, kind) a,
out size_t infoResult,
out bool lowApiExecuted,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why this flag is needed?

Copy link
Author

@aferust aferust May 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

with the original code, enforce was only triggered in the else block of if (m == 0 || n == 0), meaning that any lapack function ran. With that flag it provides the same logic. I was not going to put that flag. Then that idea seemed meaningful.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

either way we check info. Do you think the flag is redundant?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Already removed

if(lowApiExecuted)
{
enum msg ="pinv: pinv was not successful due to a convergence issue during SVD calculation.";
assert(!info, msg);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you need convergence validation for pinv function?
If so, then please do the same: two implementations of pinv, where one of them is has the info out argument.

return svdresult; //transposed
}

@safe pure @nogc nothrow SvdResult!T svdImpl
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think just svd name overload is OK.

@aferust
Copy link
Author

aferust commented May 11, 2024 via email

@aferust
Copy link
Author

aferust commented May 11, 2024

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

@9il
Copy link
Contributor

9il commented May 12, 2024

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

@aferust
Copy link
Author

aferust commented May 13, 2024

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
Copy link
Contributor

jmh530 commented May 13, 2024

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

@9il
Copy link
Contributor

9il commented May 14, 2024

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
Copy link
Contributor

jmh530 commented May 14, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants