Get-MsolUserByStrongAuthentication does not accept "Disabled" state

I want to find all Azure AD users that do not have MFA enabled
https://docs.microsoft.com/en-us/powershell/module/msonline/get-msoluserbystrongauthentication?view=azureadps-1.0 looks like what I need but something seems wrong with the -Requirements argument.
I can specify "Enabled" or "Enforced" states as expected:
PS C:\> $st = New-Object -TypeName Microsoft.Online.Administration.StrongAuthenticationRequirement
PS C:\> $st.RelyingParty = "*"
PS C:\> $st.State = "Enabled"
PS C:\> Get-MsolUserByStrongAuthentication -Requirements @($st)

UserPrincipalName                         DisplayName isLicensed
-----------------                         ----------- ----------
<EMAIL_ADDRESS>Some User   True

But when I use state "Disabled" I get an error:
PS C:\> $st.State = "Disabled"
PS C:\> Get-MsolUserByStrongAuthentication -Requirements @($st)
Get-MsolUserByStrongAuthentication : Invalid value for parameter.  Parameter name: StrongAuthenticationRequirements.
At line:1 char:1
+ Get-MsolUserByStrongAuthentication -Requirements @($st)
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    + CategoryInfo          : OperationStopped: (:) [Get-MsolUserByStrongAuthentication], MicrosoftOnlineException
    + FullyQualifiedErrorId : Microsoft.Online.Administration.Automation.PropertyValidationException,Microsoft.Online.
   Administration.Automation.GetUserByStrongAuthentication

I also tried using an empty string ""  and $null because when I get users with Get-MsolUser I notice the MFA Status is blank rather than the string "Disabled" but this value seems to give me the "Enabled" result set again.
How do I return MFA Disabled users with this cmdlet please?

@james-stocks
Have you tried:
Get-MsolUser -all | select DisplayName,UserPrincipalName,@{N="MFA Status"; E={
if( $.StrongAuthenticationRequirements.State -ne $null) {$.StrongAuthenticationRequirements.State} else { "Disabled"}}

Hi @MarileeTurscak-MSFT -
Thank you very much for looking into this!
The query you suggest works, however I need this to work on very large (>100K) user directories and I'm concerned that Get-MsolUser will cause huge amounts of data to be copied about. I was particularly interested in Get-MsolUserByStrongAuthentication because that might provide an option that only works with the data of users who do not have MFA enabled.

Please make sure you have finished the following steps 1-3 and run the following Windows PowerShell again:

Install Microsoft Online Services Sign-In Assistant for IT Professionals RTW

Install Azure Active Directory Module for Windows PowerShell (64-bit version) (64-bit is recommended)

Connect to Windows Azure AD by running the PowerShell command import-module MSOnline.

Then you can try running the commands again.

@yugangw-msft

Hi @james-stocks , are you still having this issue?

Hi Marilee-
I am still unable to select only MFA-disabled users with the Get-MsolUserByStrongAuthentication cmdlet as per the original issue.
I can certainly use Get-MsolUser and the command you provided works in my dev environment but I am unsure if running Get-MsolUser -all then piping the results to Select is efficient in Production where we have over 100,000 users in AD and only need to report on 100 non-compliant users.
I attempted the instructions in your last comment but the link for step 2 leads to a Microsoft Connect Has Been Retired page. I used the page https://docs.microsoft.com/en-us/office365/enterprise/powershell/connect-to-office-365-powershell for installation instructions, and I've run Update-Module on MSOnline and AzureAD to try to ensure I have the latest versions of these modules.

The StrongAuthenticationRequirements are applicable for "Per-User MFA" so a user in "Disabled" state can still use MFA via Conditional Access. The better option would be to look for the StrongAuthenticationMethods property which tells if the user registered any methods.
Sample: Get-MsolUser -All | where {$_.StrongAuthenticationMethods.Count -eq 0} | Select-Object -Property UserPrincipalName | Sort-Object userprincipalname
I'm not sure if there's really a more efficient option that this.

Thanks Marilee.
It's interesting that that query is users with MFA registered, and not just uses with MFA enabled. I think we need compliance checks for both.
Is there a better forum/helpdesk for requesting new interfaces for MSOL management?
Thanks!

As of today the options are limited, but you can try to scope it by user other parameters than "-All".
You can leave feature requests in office365.uservoice.com: https://office365.uservoice.com/forums/264636-general/suggestions/32706778-get-msoluser-returns-data-for-all-users
Or in the AAD feedback page: https://feedback.azure.com/forums/169401-azure-active-directory

I am closing this out for now but will update the thread upon new releases.

"Get-MsolUserByStrongAuthentication -Requirements @()" works for me


[zh-cn]sync controller-manager-leader-migration

content/zh-cn/docs/tasks/administer-cluster/controller-manager-leader-migration.md

/lgtm
/approve


FIX link to install guide

The link to the installation guide in getting_started_with_selene.ipynb is broken.

Thanks very much for catching that @VarIr!


Is it possible to specify a custom partition scheme?

I could only find this method in BrokerServiceLocator that creates partition keys for connecting to the broker service:
https://github.com/SoCreate/service-fabric-pub-sub/blob/5b92715c64c5384f10fd532d4912f4cf5c1f8409/src/SoCreate.ServiceFabric.PubSub/Helpers/BrokerServiceLocator.cs#L160
Is it possible to specify a custom partition key, rather than partitioning based on the message Type?  This way we could scale based on the message subject/topic as opposed to the message type.

It doesn't currently support specifying a custom partition scheme, but that does sound like a useful feature.  PR's are welcome if you want to implement it.

@danadesrosiers Thanks again, if we decide to use this in our current project I'll definitely submit a PR for this.

This feature is very attractive to me. I am trying to ball park the scope of work needed to implement it.

A new DataMember "PartitionKey" would need to be added to MessageWrapper to partition on
Update Broker Service to toggle between using PartitionKey and MessageType
Update Subscribe Attributes to support PartitionKey and MessageType
Other things I didn't think of?

This feature is very attractive to me. I am trying to ball park the scope of work needed to implement it.

A new DataMember "PartitionKey" would need to be added to MessageWrapper to partition on
Update Broker Service to toggle between using PartitionKey and MessageType
Update Subscribe Attributes to support PartitionKey and MessageType
Other things I didn't think of?

Can you explain what you hope to get out of this feature?  The most basic enhancement I can think of would be to update BrokerServiceLocator to change private async Task<ServicePartitionKey> GetPartitionForMessageAsync(string messageTypeName) to protected virtual.  That would allow the user to override it to specify which types to assign to which partition.  Probably registering some kind of PartitionStrategy with BrokerServiceLocator would be a better pattern.
But that would not give you the ability to partition a single message type across multiple partitions.  If you need that, I think it would be significantly more complex.

Can you explain what you hope to get out of this feature?  The most basic enhancement I can think of would be to update BrokerServiceLocator to change private async Task<ServicePartitionKey> GetPartitionForMessageAsync(string messageTypeName) to protected virtual.  That would allow the user to override it to specify which types to assign to which partition.  Probably registering some kind of PartitionStrategy with BrokerServiceLocator would be a better pattern.
But that would not give you the ability to partition a single message type across multiple partitions.  If you need that, I think it would be significantly more complex.

After talking with some coworkers I realize I was interested in the idea of "Named Queues" rather than changing the partition scheme. Sorry for the confusion please disregard the previous comment.

After talking with some coworkers I realize I was interested in the idea of "Named Queues" rather than changing the partition scheme. Sorry for the confusion please disregard the previous comment.


Replace nose with pytest

I'd be happy if someone could review this? It's very straightforward, but should still be looked at. If it takes too long, i'll just merge it.

Resolves #157


v5.0.0 - Missing Endpoints.

As of the v5.0.0 release a number of endpoints have been removed from this project. This project currently supports:

OAuth2Endpoint
ImageEndpoint

The official API at https://apidocs.imgur.com supports five endpoints:

Account
Comment
Album
Gallery
Image

Album, Comment, Gallery, and Account will be added back in this order in future releases.
They were removed as it had been some time since the project was updated and Imgur's API seems to have drifted a bit.
The rational behind this is I'd rather a new version with limited but working features than an old version with everything only partially working.

Ah, that explains why I couldn't find the Gallery endpoint. Cheers!

Since this is more than a year old, is it save to assume this project is now unmaintained?


Environment Variables keeps reading from .env

running ./gradlew bundleRelease to generate app bundles to release to google play store.
Debug mode reads correct environment variants with respective .env files. However, all variant's app bundles are being generated with only the .env file.
Root

app/build.gradle

tried adding proguard override -keep class com.dottid.BuildConfig { *; }
tried ENVFILE=.env.prod ./gradlew bundleRelease, but getting a bunch of gradle errors, ... without declaring an explicit or implicit dependency
tried changing name of .env to .env.developer but then getting a missing env error when i run command ./gradlew bundleRelease

Thank you for your time

are you using monterey?

monterey

I am using Bumblebee but get the same error. without .env file, build fails.

The specific issue from the topic starter is likely because he is building using ./gradlew bundleRelease to build all the release variants at one go instead of specifically specifying the full variant eg. ./gradlew bundleProductionRelease.
You can look at https://github.com/luggit/react-native-config/blob/master/android/dotenv.gradle in this repository to see the logic for how the env file is determined. From the logic I don't think it is possible for it to determine the correct env file to use if you are building more than one build variant at a time. In which case it will fallback on .env.
As a sanity check you can after building see android/app/build/generated/source/buildConfig/{productFlavor}/{buildType}/{package}/BuildConfig.java to check if the values are as expected.

I just updated Android studio & command line tools. Since then, I'm not receiving the the ENVFILE variable (with System.env['ENVFILE']) when dotenv.gradle is executed (Checked by printout output of System.env in dotenv.gradle).
Could there be a change in the new version of either android studio or command line tools which affects the ability to read this variable?

Just tried. Using Android Studio Chipmunk on Windows and tested with SET ENVFILE=.env.staging && react-native run-android in package.json and it works. That said I made sure that I didn't have any pre-existing gradle daemons running. There might be potential issues with a pre-existing gradle daemon and reading the current environment variables. I actually don't use the ENVFILE environment variable and just configure variants for Android and schemes for iOS.

That said I made sure that I didn't have any pre-existing gradle daemons running. There might be potential issues with a pre-existing gradle daemon and reading the current environment variables.

Good point. After retrying without pre-existing daemon running, the ENVFILE variable was read correctly.
This may seem pretty common sense in retrospect, but might still be worth adding a small note to the README for anyone who encounters the same issue.

Similar thing to what @Mr0cket reports happening here. The ENVFILE=.env.something way of setting the env file seems to have stopped working as expected.
Up till now I have had different npm scripts, each doing ENVFILE=.env.something react-native run-android, and the envfile would be correctly changed if I ran two different builds in a row.
Now the first run will pick the correct env file, and the ones following will use the same ENVFILE variable, unless the gradle daemon is stopped with ./gradlew --stop before the build.
Not sure if it can have something to do but I've recently updated my JDK version from JDK 8 to 11 (Zulu OpenJDK as specified in the latest RN docs) and also the Xcode Command Line Tools to 13.2.0.
I'm running in Mac M1
react-native: 0.67.2
react-native-config: 1.4.6
macOS: BigSur (11.6.2)

my solution is cp .env.prod .env before run a script, and .env was ignored from git

my solution is cp .env.prod .env before run a script, and .env was ignored from git

only this worked.
project.ext.envConfigFiles and ENVFILE=.env.something are ignored
"react-native-config": "^1.4.6"

./gradlew --stop

This one worked for me finally!

right above project.ext.envConfigFiles, set this:
project.ext.defaultEnvFile ="../path/to/.env.file"

Yes super annoying, seems to be an issue with m1 or montery keeping gradle deamon in memory. For fastlane added the -- stop command in the build script
gradle(
  task: "--stop"
)

Hi,
I had run into the same issue and did the following approach.
Instead of providing envConfigFiles at top level, what I did was provide default env file in the product flavor itself
    flavorDimensions "phase"
    productFlavors {
        dev {
            dimension "phase"
            project.ext.defaultEnvFile = ".env.dev"
            apply from: project(':react-native-config').projectDir.getPath() + "/dotenv.gradle"
            applicationId project.env.get("APP_ID")
        }
        stage {
            dimension "phase"
            project.ext.defaultEnvFile = ".env.stage"
            apply from: project(':react-native-config').projectDir.getPath() + "/dotenv.gradle"
            applicationId project.env.get("APP_ID")
        }
        prod {
            dimension "phase"
            project.ext.defaultEnvFile = ".env.prod"
            apply from: project(':react-native-config').projectDir.getPath() + "/dotenv.gradle"
            applicationId project.env.get("APP_ID")
        }
    }

This should allow to build multiple variants with one command.
Since apply is getting invoked multiple times there might be a performance hit. Maybe someone with gradle experience can comment here
One drawback here is that env variable should be used only within flavor definition. Once it goes outside the block, it may be overwritten by the next flavor.
Also more complex combination of flavors may become difficult to handle.

I just stumbled on this as well... the strange this is that it's only affecting one of my projects, not both. ./gradlew --stop worked for me so I've added it to my fastlane as per #issuecomment-1230403293

right above project.ext.envConfigFiles, set this:
project.ext.defaultEnvFile ="../path/to/.env.file"

This one line worked like a charm in my case

When running it locally the java process also stays running forever on the m1 new mac books, taking quite some memory, so works better to stop the process.
gradle(
task: "--stop"
)

from your root folder try this
cd android && ./gradlew clean && ENVFILE=../.env.production ./gradlew bundleRelease && cd ..

my solution is cp .env.prod .env before run a script, and .env was ignored from git

I have a few apps using the same codebase and was attempting to keep the .env files in separate folders. Using the ENVFILE=app1/.env.development command would not work for me at all. Only the above has worked correctly

my solution is cp .env.prod .env before run a script, and .env was ignored from git

only this worked. project.ext.envConfigFiles and ENVFILE=.env.something are ignored
"react-native-config": "^1.4.6"
"./gradlew --stop" - this helped to reset to .env from .env.beta (no idea how it .env.beta was selected before 🤨)

Yes, worked for me

Hi , I am still facing this issue !
The release builds always picks up .env config instead of any other extenstion.

add defaultEnvFile=.env.prod in to file gradle.properties

In my case this issue was caused by calling
apply from: project(':react-native-config').projectDir.getPath() + "/dotenv.gradle"
too early. It must be called AFTER:
project.ext.envConfigFiles = [     debug: ".env.development",     release: ".env.production",     anothercustombuild: ".env", ]


chore: add private ip restriction note to LangChain on VertexAI RAG tutorials

Description
Thank you for opening a Pull Request!
Before submitting your PR, there are a few things you can do to make sure it goes smoothly:

[x] Follow the CONTRIBUTING Guide.
[x] You are listed as the author in your notebook or README file.

[x] Your account is listed in CODEOWNERS for the file(s).

[x] Make your Pull Request title in the https://www.conventionalcommits.org/ specification.
[x] Ensure the tests and linter pass (Run nox -s format from the repository root to format).
[x] Appropriate docs were updated (if necessary)

Fixes #<issue_number_goes_here> 🦕

Linting error: gemini/reasoning-engine/tutorial_cloud_sql_pg_rag_agent.ipynb:cell_16:6: error: Name "IAM_USER" is not defined  [name-defined]
IAM_USER is defined. Is there a way to ignore this error? Thanks!


merge package from arlac77/npm-package-template

package.json

chore(package<EMAIL_ADDRESS>

Coverage remained the same at 80.878% when pulling 42f9bb4d539d81486ac6dffd3932e443ddb3fb14 on template-sync-1 into e0c316e1fad9d0ea9ced6623aa82c38b270a3cfd on master.

:tada: This PR is included in version 7.8.37 :tada:
The release is available on:

npm package (@latest dist-tag)
GitHub release

Your semantic-release bot :package::rocket:


Update reflection.md - acronym clarification

clarification on IL acronym, assumes a level of knowledge beyond the target audience

Thanks again @Skuldo
I've approved the changes, and I'll :shipit: now.


Feat/header icons and text

• replaced icons with format of column so you can see difference of number and date columns
• added icons for primary_key and foreign_key so it's easy to see those columns

updating green to brand color

@phamhieu i'll leave in the icon switch function, in case we have a last minute change of mind.


refactor: make every get_value function delegate to get_details

This PR
Makes sure that each get_<type>_value method in OpenFeatureClient delegates to the corresponding get_<type>_details method. Some of these methods already did so but others called evaluate_flag_details directly.
To adhere to a principle of least surprise, remaining methods were adapted to follow a common calling convention.

I'll rerun the GitHub action.


Unable to install build-essential

My Dockerfile:
FROM blitznote/debootstrap-amd64:16.04
RUN apt-get update
&& apt-get upgrade -y
&& apt-get install -y
build-essential
&& apt-get clean
&& rm -rf /var/lib/apt/lists/*
The Error:
The following packages have unmet dependencies:
build-essential : Depends: libc6-dev but it is not going to be installed or
libc-dev
Depends: g++ (>= 4:5.2) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

Confirmed, can reproduce. Images in staging have already moved past this.

I've edited your report, minimized the supplied dockerfile and added formatting.
BTW, don't use apt-get upgrade in dockerfiles. It's a bad practice.
Any derived image should be rebuild on rebuilds of its parent image(s).

Thanks.  I don't normally use upgrade.  I just added it to see if it made a difference for this issue and left it there, my bad.

New baseimages are up.


Adjust http client test timeouts

Motivation
Tests are failing sporadically on Travis for the operations maintainer. So increase the timeouts that they should not fail because of timeouts

Coverage decreased (-0.07%) to 86.816% when pulling 92742e6523bb3580da83ac63f6547f1a4a15727e on test_timeouts into 8a1cc3932b41629a7c455a8daec2ad5d508c7e53 on master.


bug(Configuration.java) configuration builder gets defaults when not needed

Configuration builder forces a dependency on json-smart when defaults are gotten, even when they are not needed.

fixes #121

Fixed by #128


Fixed recalculation after lane connection remove

fixed recalculation after lane connection removal,
refactored RequestRecalculation(segmentId, propagate)

Closes #1198

its working!


NEXT-00000 - Return 401 status code on failed login

1. Why is this change necessary?
A failed login in the storefront will return a response with the statuscode 200. This should be a 401.
A failed login in the administration will return a response with the statuscode 400. This should be a 401.
2. What does this change do, exactly?
This changes the statuscode for a failed login attempt in the storefront/administration to a 401.
3. Describe each step to reproduce the issue or behaviour.
Go to the default administration login. Fail your login and check the statuscode. It is a 400.
Go to the default storefront login. Fail your login and check the statuscode. Its is a 200.
4. Please link to the relevant issues (if any).
/
5. Checklist

[x] I have rebased my changes to remove merge conflicts
[ ] I have written tests and verified that they fail without my change
[x] I have created a changelog file with all necessary information about my changes
[ ] I have written or adjusted the documentation according to my changes
[ ] This change has comments for package types, values, functions, and non-obvious lines of code
[x] I have read the contribution requirements and fulfil them.

@jasperP98 are you still willing to work on this? :slightly_smiling_face:

@mitelg I agree with the approach above but do not have the time right now to implement such change

@jasperP98 no worries :slightly_smiling_face: :v:  to keep the PRs clean, I will close this one for now, and if you find the time again, feel free to open up a new PR with the comments considered :+1:


Dynamic proxy

Dynamic Proxy.
Supports Etcd, Consul and ZooKeeper using github.com/docker/libkv.
Demonstration Video: https://youtu.be/I0Kax0F1XWM
Examples
Etcd
proxy / etcd://<IP_ADDRESS>:2379/caddyserver.com

Consul
proxy / consul://<IP_ADDRESS>:8500/caddyserver.com

So after watching your video, this definitely does seem to be shaping up! I'm really glad you've taken the plunge into getting support for dynamic backends at least started!
Based on some experience and seeing what others provide, I have some input that may not exactly apply to this specific PR but can at least spawn discussion and thinking for future additions to this end goal.

It would be really cool if Caddy could provide the ability to consume proxy based configuration via container runtime labels or variables. You could consume a label/variable like caddy.port=2015 telling the proxy to utilize the public port entry in etcd mapped to 2015. You could even define http or https backend connectivity i.e caddy.protocol=https or weight caddy.weight=7 for this particular backend.

Note that this depends on what registrator dumps in your K/V store based on metadata it finds from your docker socket! I know for a fact a lot of things are possible with environment variables.

Not sure much about the etcd entry with the appended /caddy-2015 but to me, the only data in the actual front-end proxy configuration (relating to the container infra) should be the name of the service. In example, if I have a service bugzilla, I should only really need to reference /bugzilla in etcd. That's, for all intensive purposes, the key for which back-end containers register themselves with a front-end.

Global definition of backend store locations is definitely something I'd like to see. A lot of the time, docker infrastructure has multiple redundant proxies all making backend connections to a multitude of container hosts. For an entire deployment like this, the proxies would essentially utilize the same K/V store and service discovery system for the entire cluster. Maybe utilize a variable or something for this i.e. proxy / main_etcd_store://bugzilla with something like dynamic etc etcd://<IP_ADDRESS>:4444 setup globally. Keep reading for more ideas in this realm.

Re-iteration on the last point, why don't you allow for sub-service declarations under a specific from and even specify the destination K/V store right in the stanza i.e.

proxy / {
  bugzilla {
    destination: etcd
    max_fails 7
    health_check /bugs
  }
  jira {
    fail_timeout 30
    destination: consul
  }
}

So this says look in the globally defined etcd K/V for bugzilla and health check any backend container endpoints @ /bugs to ensure a proper response is returned. Likewise for jira, just a tad different.

I haven't used Caddy extensively yet but it seems it doesn't yet support Name Based Proxying? This feature is a must for Caddy to compete long-term as a docker load-balancer and reverse proxy.
Something like this would be very exciting:

proxy / {
  bitbucket {
    name: git.example.com
    destination: zookeeper

I'm subscribed to this thread so post replies and lets work together to figure out the best way to do this!

....and now I am subscribed as well.  Is there an IRC or messaging client that the three of us have in common?

@faddat Join us on gitter for dev chat: https://gitter.im/mholt/caddy

Alright, having a go at this now.  I'll send feedback when I'm finished up with it.  I spoke with @abiosoft not long ago with the intent of doing it then, but was thwarted-- off to server-land with me!

What does this need to be mergable?

After Caddy 0.9 is released (not next week, but the release after -- not sure when that is yet), Caddy will be in a much better position to have this merged. Not quite sure of the details yet but I'll work closely with @abiosoft at that point.

@mholt What will be the performance impact of this PR? Benchmark?

@thalesfsp No clue, honestly - @abiosoft might know better. I would imagine it would only have minimal, if any, impact, since as I understand it, the updating of proxies is done in the background.

@mholt Awesome! I hope and will wait for @abiosoft come with some benchmark :)

You can always perform your own, too.

I expect it to have little to no impact

It's amazing this PR can still be merged :smile: That will not be the case for long though.
The 0.9 changes are in a branch and will soon be in master. I don't yet have the proxy middleware brought over to that branch but I will soon, and then these changes can possibly go on top.

I use caddy as reverse proxy with docker-gen for my docker host and it works great. But etcd support sounds like a better solution. Can't wait to see this feature and the howto to use it ;)

@mholt  Whats the plan to when merge this PR?

I think it depends on @abiosoft's schedule.

Would this lay the groundwork for caddy being able to supply an ingress service to kubernetes? As far as I see it... an ingress service is basically a dynamic proxy using the k8s apiserver instead of general key-value stores.

@abiosoft @mholt What's the ETA to be merged?

No ETAs, sorry. I plan to resume work on this after 0.9 :)
There is a very good chance it will make it in before 1.0 and could even be earlier depending on Caddy's priorities.

A+
Would really help with Rancher / Kubernetes environments.

Yes, would use it with rancher...

@abiosoft Do you have an updated branch? I have interest in use it. :)

@abiosoft 0.9 released ;)

@mholt @abiosoft 0.9 released ;)

Okay, good news :)
Work is about to resume on this. However, I wanna make this a pluggable feature so proxy will only be making use of it, rather than having it built into proxy. This will make it easy to bring the feature to any other directive or plugin that needs it.
Still, no ETA 😄 but I will keep this thread updated.

@abiosoft great idea! much better and extensible way to implement.
Have you had anymore thoughts on implementation details? I'm wondering if it's not worth getting this working the original intended way (when you created this PR), and then move on to other ways later (since this PR is quite large both in conversation and commits).
Additionally, I've had some other ideas about how this could work based on other projects that have come along (i.e. Fabio, Traefik, etc..) and attempted to provide LE functionality with dynamic backends.
I would love to talk through some of this directly with you guys @abiosoft @faddat @stp-ip and others..
Seems to be a #caddy on freenode that I just joined. Ping me and we can chat (notifications go to my phone). Or maybe just create a new forum thread; whatever's easier.

@InAnimaTe I'm rarely on irc :( What about Gopher's slack ?

@abiosoft That works! I just joined. Go ahead and make a channel and invite me plz ;)

Is this PR still going? Also, will this support native docker swarm?

If this supported swarm mode,  wow.   Also,  I think that it isn't too
tough to implement.   Right now I haven't time,  but it's possible that
myself or someone on my team will do it.
... Because Caddy is awesome.
On Dec 7, 2016 9:25 PM, "HazCod"<EMAIL_ADDRESS>wrote:

Is this PR still going? Also, will this support native docker swarm?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/mholt/caddy/pull/564#issuecomment-265459980, or mute
the thread
https://github.com/notifications/unsubscribe-auth/AGz6iXhD7LaOTu2KYRp7wEQR95ePV4uaks5rFsHggaJpZM4HRNBY
.

Maybe it is time to write it once again? :)

I've put some thought in this and isn't it cleaner to do this via docker-gen ?
Altough I am unsure how it works for swarm/multiple hosts.

I think I'm going to close this PR -- it's almost a year old, but I hope it will soon make it to Caddy -- if not into core, then maybe it would do well as a plugin? Thanks for all the work you've put into it, @abiosoft, I hope we can use it someday. :)


A different approach to compute besselk in intermediate zones

@cgeoga I was thinking some more about how we compute the region between asymptotic expansions (power series, uniform expansions, and large argument). Which corresponds to roughly 0<nu<25 and 2 < x < 30...
Right now, we are using the wronskian relation to besseli. The reason we can do this is that the power series of besseli is convergent in this zone so using that along with the continued fraction and wronskian we can get besselk pretty accurately. Now there are a couple drawbacks to that:

We need besseli(nu, x) and besseli(nu + 1, x) and then also need to compute besselk(nu+1,x)/besselk(nu,x).
Power series are slower to converge as x gets bigger
The power series of besseli(nu, x) suffers from cancellation for complex numbers so can't be as uniformly used.

The first point was taken care of by using SIMD instructions to compute both besseli(nu, x) and besseli(nu + 1, x). So that is pretty optimized (though we pay an extra penalty for complex numbers). But points 2 and 3 are pretty much unsolvable. Power series will always be slow for larger arguments and can't use much for complex numbers.
I was toying with an alternative which links a few ideas. Now the continued fraction approach to compute besselk(nu+1,x)/besselk(nu,x)converges more quickly asxgets bigger so it has an inverse relation to the power series. (i.e., power series converges faster for smallxwhile continued fraction converges faster for bigx`). So what if there was a good way to use the uniform expansions to link with the smaller orders. Now Miller's algorithm links stable backward recurrence to higher orders with a normalization factor calculated from low orders. That is fine but usually we do not know the lower orders unless for integer arguments.
I'm sure this idea has been described before but can't seem to find a description of it. But there appears to be a way to link stable forward recurrence with a normalization factor calculated at higher orders to compute besselk at lower orders. It works by instead of considering trial start values (0 and epsilon) in the Miller backwards scheme, it uses the continued fraction approach to generate start values in the forward scheme. So the algorithm would look something like:

Generate start values 1.0 and knup1/knu.
Use upward recurrence up until a order you can use asymptotic expansion (nu_ceil).
Compute besselk(nu_ceil, x) using asymptotic expansion
Normalize starting value

It may seem complicated but instead of having to compute besseli(nu, x), besseli(nu + 1, x) and besselk(nu+1,x)/besselk(nu,x)we just needbesselk(nu_ceil, x)andbesselk(nu+1,x)/besselk(nu,x). There is a few extra steps of course to generate the recurrence but that is pretty small.
So now we have an algorithm that computes besselk(nu_ceil, x) at constant time using uniform asymptotic expansion and algorithm to compute besselk(nu+1,x)/besselk(nu,x)which converges faster asx` increases. A sample implementation of the algorithm would look like..
using Bessels

function test(nu, x)
    fpart, _ = modf(nu)
    nu_ceil = 25.0 + fpart
    knu_ceil = Bessels.besselk_large_orders(nu_ceil, x)
    knup1 = Bessels.besselk_ratio_knu_knup1(nu, x)
    return knu_ceil / Bessels.besselk_up_recurrence(x, knup1, 1.0, nu + 1, nu_ceil)[1]
end

Speed tests to other approach.
## high precision calculation
julia> 2.308550810254649099632E-14
2.308550810254649e-14

julia> test(2.2, 30.0)
2.3085508102546522e-14

julia> Bessels.besselk_continued_fraction(2.2, 30.0)
2.3085508102546462e-14

julia> @benchmark test(2.2, x) setup = (x = 30.0 + rand())
BenchmarkTools.Trial: 10000 samples with 968 evaluations.
 Range (min … max):  79.846 ns … 120.610 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     79.976 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   80.832 ns ±   2.255 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆█▃                    ▂▃           ▃▅▁                      ▁
  ███▆▅▅▄▃▄▄▄▄▆▄▄▃▁▁▄▃▄▁▆██▆▅▃▃▃▄▄▄▅▄▆███▅▄▅▃▃▁▃▁▄▃▁▄▃▃▁▃▃▁▁▃█ █
  79.8 ns       Histogram: log(frequency) by time      87.2 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 30.0 + rand())
BenchmarkTools.Trial: 10000 samples with 898 evaluations.
 Range (min … max):  127.180 ns … 178.360 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     129.315 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   128.840 ns ±   1.969 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▇▂              █▇
  ▃███▃▂▂▂▂▂▂▂▂▁▂▂▂███▄▃▂▂▂▂▂▃▃▃▄▃▂▂▂▂▂▂▁▂▂▂▂▃▃▂▂▂▂▂▁▁▂▂▂▂▂▂▂▂▂ ▃
  127 ns           Histogram: frequency by time          134 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Which is about 40% faster for bigger arguments. For smallish arguments its about 10% faster.
julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 10.0 + rand())
BenchmarkTools.Trial: 10000 samples with 940 evaluations.
 Range (min … max):  102.305 ns … 145.612 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     104.122 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   104.488 ns ±   2.483 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅▆     ▄█▅▂   ▁    ▅▄▃▁     ▁                                 ▂
  ██▇▁▁▁▁████▆▄▄█▇▅▅▄████▇▅▅▄▇█▇▆▆▇▇▅▄▁▃▄▅▅▆▅▆▅▆▆▅▄▅▄▄▃▄▄▅▅▃▄▄▄ █
  102 ns        Histogram: log(frequency) by time        115 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark test(2.2, x) setup = (x = 10.0 + rand())
BenchmarkTools.Trial: 10000 samples with 954 evaluations.
 Range (min … max):  93.203 ns … 134.565 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     97.615 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   96.402 ns ±   2.687 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇▅              ▁       █▅              ▁         ▁▃         ▂
  ██▇▅▄▅▄▅▅▄▄▁▁▃▄▅█▆▅▄▁▁▁███▇▅▅▅▅▇▆▆▅▁▄▁▃▇█▇▆▅▅▄▃▅▁▅██▆▃▄▁▄▁▄▃ █
  93.2 ns       Histogram: log(frequency) by time       104 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And for much smaller arguments it's hard to beat the power series so it is slower.
julia> @benchmark test(2.2, x) setup = (x = 4.0 + rand())
BenchmarkTools.Trial: 10000 samples with 918 evaluations.
 Range (min … max):  115.423 ns … 157.588 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     121.642 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   122.518 ns ±   4.119 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆              █             ▂█              ▆              ▃ ▂
  █▇▃▄▄▆▃▃▁▁▃▁▁▁▇██▅▃▆█▆▄▄▄▃▃▄▁███▄▁▅█▆▄▃▄▁▃▃▄▅██▃▄▄▅▅▃▁▃▁▄▁▃▁█ █
  115 ns        Histogram: log(frequency) by time        133 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 4.0 + rand())
BenchmarkTools.Trial: 10000 samples with 942 evaluations.
 Range (min … max):  100.849 ns … 158.882 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     108.678 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   108.736 ns ±   4.641 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄    ▅      █▂  ▅▂  ▂▆  ▇▄   ▃  ▂▅  ▅▃   ▁   ▃   ▃            ▂
  █▄▁▄▇█▅▁▆▃▃▄██▆▆██▄▄███▅██▆▆██▆▄██▇▆██▅▆▅█▆▅██▇▅▆█▅▅▄▇▄▄▄██▁█ █
  101 ns        Histogram: log(frequency) by time        123 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

So there is of course things to consider but this approach has the advantage of being more applicable to a large parameter set including complex numbers and has the potential for large speed gains over a fairly big domain. My initial implementation could also probably be improved to speed things up a bit but it does serve as a good way to link directly to the asymptotic expansions.

Not a problem!! I was working on this before I needed to switch to dissertation stuff so haven't come back to it. Hopefully will work on it a little more in the next couple of weeks.
Yes - the power series for besseli is pretty optimized and rapidly convergent for small arguments so it is hard to beat for x<3. But I was also thinking about the idea of using this approach for bessely because you can't use the besselj power series the same way you can for besseli. So it is is much more applicable to that function but also here when the besseli is slowly converging (larger arguments) or prone to severe cancellation (complex numbers).
I also like the idea of reducing the number of algorithms employed such that if the uniform expansion is the main driver of the speed and accuracy of almost all arguments that is more easily optimized than several different routines. But ya, the Miller algorithm seems to be a pretty common way to generate higher order (integer orders) functions when forward recurrence is not stable. That way you can use an optimized zero order function to normalize for any order. This doesn't work as well for besselk and bessely though when forward recurrence is stable..
Anyway, still need to have the continued fraction approach for bessely..... that could reduce significantly the amount of code and compilation costs (though that seems to be less of an issue now)

Was thinking about this approach for bessely again.... it might be possible to derive the continued fraction for bessely but I think its also possible to use the relation between besselk and the Hankel function. Then use the relation of the Hankel function to get bessely from besselk therefore we can just use the same approach described above using besselk_ratio_knu_knup1 with some modifications.
The first challenge is that it requires a rotation of argument by pi/2 so real arguments for bessely will require evaluation of besselk along imaginary axis. That's not too big of a deal for the continued fraction approach other than slowing it down slightly. I made some adjustments to the existing approach here to speed it up slightly for complex arguments....
function besselk_ratio_knu_knup1(v, x::ComplexOrReal{T}) where T
    MaxIter = 1000
    S = eltype(x)
    (hn, Dn, Cn) = (S(1e-50), zero(S), S(1e-50))

    jf = one(T)
    vv = v * v
    @fastmath begin
        for _ in 1:MaxIter
            an = vv - (jf - 0.5)^2
            bn = 2 * (x + jf)
            Cn = an / Cn + bn
            Dn = inv(muladd(an, Dn, bn))
            del = Dn * Cn
            hn *= del
            checktol(del - 1) && break
            jf += one(T)
        end
        xinv = (v + x + hn + 1/2) / x
    end
    return xinv
end
@inline checktol(err::T) where T = abs(err) < eps(Float64)
@inline checktol(err::Complex{T}) where T = abs2(err) < eps(T)^2

Perhaps, it could be sped up further but it is pretty fast already for the input range we need for bessely (7<x<20)...
julia> @btime besselk_ratio_knu_knup1(v, x) setup=(x= 20.0im; v = rand()*1)
  76.246 ns (0 allocations: 0 bytes)
1.0000915788150249 - 0.045999085366594im

julia> @btime besselk_ratio_knu_knup1(v, x) setup=(x= 7.0im; v = rand()*1)
  141.996 ns (0 allocations: 0 bytes)
1.001602364037009 - 0.11366618431295909im

So a range between 76-142 nanoseconds in the regime we'd utilize this for even with complex arithmetic... Of course we also need recurrence up to higher order expansion (takes around 20 nanoseconds) and then need one large order evaluation...
julia> @btime Bessels.besseljy_debye(y, x) setup=(x=10*rand()+12.0; y=rand()+60.0)
  49.393 ns (0 allocations: 0 bytes)
(1.8133319226273097e-26, -3.0253831328694315e23)

Which takes around 50 nanoseconds. So that puts the total time for this approach around 150-225 nanoseconds in a region where no expansion gives good results. Now that might not seem that fast but the existing approach using Chebyshev approximation is the longest approach we have...
julia> @btime Bessels.bessely_chebyshev(y, x) setup=(x=12*rand()+6; y=rand()*20.0)
  460.577 ns (0 allocations: 0 bytes)

Which takes between 460-500 nanoseconds so the new approach would be about 2-3x faster.
I think it might be worth taking a look at the whole routine to some degree and try to decrease the amount of algorithms to improve branch prediction.
However, there was some discussion on discourse about this but it might be possible to improve the Chebyshev evaluation using SIMD.... I'd assume it would be possible to improve those times by 2-3x

@cgeoga To continue this discussion a little bit because I think it would be really nice to round out the modified bessel functions completely. I've also been thinking a lot about how to ensure the implementations are differentiable. Enzyme has come a long way on the main branch and can give pretty fast derivatives of the code in Bessels.jl already.
Anyway for computation of modified bessel functions (real arguments for now) there are really three main approaches in order of precedence.

large argument expansion
power series
uniform large argument expansion for large orders

Now the power series for besseli is always numerically stable. And with some of the recent work with SIMD we should probably unroll this loop by a factor of four and use SIMDMath.jl to make sure that vectorizes. That will make computation of the power series much more efficient. But where all of those methods fail (let's call it the intermediate zone) we need a good fourth algorithm. At first it doesn't seem like there is a great method but this thread mentions two:

besseli power series coupled with continued fraction for besselk(v+1)/besselk(v) and wronskian.
The slightly modified Miller's scheme to use with forward recurrence and continued fraction with normalizing condition with large argument expansion

Those methods advantages and disadvantages were extensively discussed above but essentially method 1 converges faster for smaller arguments while method two converges faster for larger arguments. The slight disadvantage of method 2 is that the accuracy is tied to the uniform expansions which is slightly less accurate (due to roundoff) than the power series. Though approach 1 is not feasible for complex numbers.
I'm going to introduce a different approach that could be a more robust approach meaning that it will work for complex numbers and be accurate. This was inspired by https://arxiv.org/pdf/1112.0072.pdf when I was working on the airy functions but approach three is:

Sequence transformation (Levin's approach) of the asymptotic expansions

At first I ignored these approaches because they honestly just didn't seem worth it and not efficient at all. Essentially they are able to sum a diverging asymptotic series and increase the validity of the asymptotic expansions. Computationally they can be computed using recursive schemes but they require significant more computation. Essentially you must compute a Levin table and compute down the counter diagonals of the Levin table. Now just thinking theoretically about this problem this requires at a minimum N*(N+1)/2 muladd operations. I ran some tests and typically you need about 16-24 terms to get this to work. So for 16 terms you need 136 muladd operations and for 24 terms you need 300 muladds. And this is at a minimum. Not to mention the sequence transformation is essentially a more advanced Pade approximant so you must also compute a numerator and denominator so all those operations need to be doubled.
This is probably why I initially avoided this approach was that most of routines are under 100 ns and not only do you need all these arithmetics you also need to store tables and so your cache lines and memory accesses are honestly rubbish my first attempt at this using arrays took like 10 microseconds which is awful. But I want to come back to the performance after discussing where this is accurate.
So here are a few plots of the relative errors as a function of x at a specific nu value. I'll then plot the number N of terms used in each sequence transformation.

So the really cool thing here is we are summing the asymptotic expansion for large x but this sequence transformation allows us to sum a rapidly diverging series even when x << nu.
Ok for the benchmarks it's hard to say as the general form was so slow but I completely unrolled a version (which is probably not what we want) for N=16 but to really see the limits of the performance. Here is the code it is long (@oscardssmith maybe if you have an idea on how to organize the memory most effectively in a loop here 😄 ) ill work on optimizing the general case in a bit. But the benchmarks are pretty good.
julia> p = ntuple(i -> i*2.0, 18)
(2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0)

julia> @benchmark g1($p)
BenchmarkTools.Trial: 10000 samples with 990 evaluations.
 Range (min … max):  45.328 ns … 84.385 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     45.455 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   45.635 ns ±  1.419 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃█▆▄▂                                   ▁                   ▁
  █████▄▁▁▁▁▁▁▃▁▁▁▃▁▄▁▃▃▁▃▅▃▁▁▁▁▃▃▁▄▁▁▁▁▃▇██▇▄▄▁▃▃▁▃▃▄▁▁▄▃▃▅▅ █
  45.3 ns      Histogram: log(frequency) by time      49.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And the advantage is this isn't awful for complex numbers because we have only real/complex multiplies and complex/complex additions.
julia> p = ntuple(i -> i*2.0 + rand()*im, 18)
(2.0 + 0.20671806918378222im, 4.0 + 0.5544791253581327im, 6.0 + 0.8743485225845774im, 8.0 + 0.5842694164427652im, 10.0 + 0.5376242211139655im, 12.0 + 0.24931450428480173im, 14.0 + 0.7889137184169184im, 16.0 + 0.18129161121223591im, 18.0 + 0.9320649223876125im, 20.0 + 0.4467818886051006im, 22.0 + 0.44671130589145225im, 24.0 + 0.09229756644375375im, 26.0 + 0.09649719037936122im, 28.0 + 0.7052070903007531im, 30.0 + 0.8544141179769782im, 32.0 + 0.35520748181052775im, 34.0 + 0.22437122283477695im, 36.0 + 0.9058299917597746im)

julia> @benchmark g1($p)
BenchmarkTools.Trial: 10000 samples with 832 evaluations.
 Range (min … max):  147.886 ns … 208.534 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     148.689 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   149.701 ns ±   3.756 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▇█▆▂     ▁▃▄▂▁▁ ▁▁                                           ▂
  █████▇▇█▇▇██████████▇▆▅▆▆▅▆▆▇▆▅▅▅▃▆▆▆▄▅▄▅▅▄▅▁▅▆▅▃▅▆▆▆▄▆▅▅▅▆▄▅ █
  148 ns        Histogram: log(frequency) by time        168 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

const CC =  ((1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
(0.5, 0.6, 0.6666666666666666, 0.7142857142857143, 0.75, 0.7777777777777778, 0.8, 0.8181818181818182, 0.8333333333333334, 0.8461538461538461, 0.8571428571428571, 0.8666666666666667, 0.875, 0.8823529411764706, 0.8888888888888888),
(0.32, 0.4166666666666667, 0.4897959183673469, 0.546875, 0.5925925925925926, 0.63, 0.6611570247933884, 0.6875, 0.7100591715976331, 0.7295918367346939, 0.7466666666666667, 0.76171875, 0.7750865051903114, 0.7870370370370371),
(0.23148148148148148, 0.31486880466472306, 0.3828125, 0.438957475994513, 0.486, 0.5259203606311045, 0.5601851851851852, 0.589895311788803, 0.6158892128279884, 0.6388148148148148, 0.6591796875, 0.6773865255444739, 0.693758573388203),
(0.17992503123698458, 0.251220703125, 0.31214753848498705, 0.3645, 0.4098080732190424, 0.44931520061728397, 0.48401666608312033, 0.514707413577676, 0.5420246913580247, 0.5664825439453125, 0.5884987009255157, 0.6084152568206066),
(0.14654541015625, 0.20809835898999135, 0.26244, 0.3104606615295776, 0.35303337191358025, 0.3909365379902126, 0.4248378651752246, 0.45530074074074073, 0.48279762268066406, 0.5077243694259351, 0.5304133008179648),
(0.12331754606814303, 0.177147, 0.22578957202151098, 0.2696782702117627, 0.30931242566258577, 0.34518076545487, 0.3777309849108368, 0.4073604941368103, 0.4344165727708536, 0.4592004039488861),
(0.1062882, 0.15394743546921202, 0.1977640648219593, 0.237932635125066, 0.27473571128040675, 0.30848030434385004, 0.3394670784473419, 0.3679763910529583, 0.39426297308742747),
(0.09330147604194668, 0.13596279456509702, 0.17570409978466411, 0.21259310991936234, 0.24678424347508002, 0.27846908778883517, 0.30784952976979524, 0.33512352712431337),
(0.08308837445644818, 0.12164129985092131, 0.15792631022581202, 0.19194330048061778, 0.22376980268745683, 0.25352314216336075, 0.28133826968460873),
(0.07485618452364388, 0.10998439462154766, 0.14331766435886129, 0.17482015834957565, 0.20452287098892968, 0.2324948200865864),
(0.06808557762286284, 0.1003223650512029, 0.13111511876218174, 0.1604100948932782, 0.18821009245104614),
(0.06242280492074847, 0.09219031787965903, 0.12077936556670357, 0.14812831350313815),
(0.057618948674786896, 0.08525602275296724, 0.11191917020237105),
(0.053493975060685324, 0.07927607889334616),
(0.049914568192106844,))

function g1(s)
    b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12, b13, b14, b15, b16, b17 = 1 / (s[2] - s[1]), 1 / (s[3] - s[2]), 1 / (s[4] - s[3]), 1 / (s[5] - s[4]), 1 / (s[6] - s[5]), 1 / (s[7] - s[6]), 1 / (s[8] - s[7]), 1 / (s[9] - s[8]), 1 / (s[10] - s[9]), 1 / (s[11] - s[10]), 1 / (s[12] - s[11]), 1 / (s[13] - s[12]), 1 / (s[14] - s[13]), 1 / (s[15] - s[14]), 1 / (s[16] - s[15]), 1 / (s[17] - s[16]), 1 / (s[18] - s[17])
    a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15, a16, a17 = s[1] * b1, s[2] * b2, s[3] * b3, s[4] * b4, s[5] * b5, s[6] * b6, s[7] * b7, s[8] * b8, s[9] * b9, s[10] * b10, s[11] * b11, s[12] * b12, s[13] * b13, s[14] * b14, s[15] * b15, s[16] * b16, s[17] * b17

    index = 1

    a1 = muladd(a1, -CC[index][1], a2)
    b1 = muladd(b1, -CC[index][1], b2)

    a2 = muladd(a2, -CC[index][2], a3)
    b2 = muladd(b2, -CC[index][2], b3)

    a3 = muladd(a3, -CC[index][3], a4)
    b3 = muladd(b3, -CC[index][3], b4)

    a4 = muladd(a4, -CC[index][4], a5)
    b4 = muladd(b4, -CC[index][4], b5)

    a5 = muladd(a5, -CC[index][5], a6)
    b5 = muladd(b5, -CC[index][5], b6)

    a6 = muladd(a6, -CC[index][6], a7)
    b6 = muladd(b6, -CC[index][6], b7)

    a7 = muladd(a7, -CC[index][7], a8)
    b7 = muladd(b7, -CC[index][7], b8)

    a8 = muladd(a8, -CC[index][8], a9)
    b8 = muladd(b8, -CC[index][8], b9)

    a9 = muladd(a9, -CC[index][9], a10)
    b9 = muladd(b9, -CC[index][9], b10)

    a10 = muladd(a10, -CC[index][10], a11)
    b10 = muladd(b10, -CC[index][10], b11)

    a11 = muladd(a11, -CC[index][11], a12)
    b11 = muladd(b11, -CC[index][11], b12)

    a12 = muladd(a12, -CC[index][12], a13)
    b12 = muladd(b12, -CC[index][12], b13)

    a13 = muladd(a13, -CC[index][13], a14)
    b13 = muladd(b13, -CC[index][13], b14)

    a14 = muladd(a14, -CC[index][14], a15)
    b14 = muladd(b14, -CC[index][14], b15)

    a15 = muladd(a15, -CC[index][15], a16)
    b15 = muladd(b15, -CC[index][15], b16)

    a16 = muladd(a16, -CC[index][16], a17)
    b16 = muladd(b16, -CC[index][16], b17)

    index = 2

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]

    a14 = a15 - a14*CC[index][14]
    b14 = b15 - b14*CC[index][14]

    a15 = a16 - a15*CC[index][15]
    b15 = b16 - b15*CC[index][15]

    index = 3

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]

    a14 = a15 - a14*CC[index][14]
    b14 = b15 - b14*CC[index][14]

    index = 4

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]

    index = 5

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    index = 6

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    index = 7

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    index = 8

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    index = 9

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    index = 10

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    index = 11

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    index = 12

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    index = 13

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    index = 14

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    index = 15

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    index = 16

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    return a1 / b1

end

Isn't this just 16 pairs of evalpolys of degrees 1 to 16?

Yeah, my suggestion was the hybrid of this looped verison with your unrolled version (where you precompute cst but otherwise leave it looped). I think that should vectorize well and be fast.

Ya I think my initial attempt at this removing that cst steered me away from this approach considering the n(n+1)/2 growth in operations for the sum. Just an example code..

function g(s)
    n = length(s) - 1
    l = similar(s)
    N = similar(l)
    D = similar(N)

    for k in 1:n
        w = s[k+1] - s[k]
        N[k], D[k] = s[k] / w, 1.0 / w
        for j in 1:k-1
            N[k-j] = muladd(-2.1,N[k-j],N[k-j+1])
            D[k-j] = muladd(-2.1,D[k-j],D[k-j+1])
        end
        l[k] = N[1]/D[1]
    end
    return l
end

julia> p = rand(16);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 200 evaluations.
 Range (min … max):  406.665 ns …  2.451 μs  ┊ GC (min … max): 0.00% … 82.54%
 Time  (median):     416.250 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   423.480 ns ± 92.645 ns  ┊ GC (mean ± σ):  1.09% ±  4.06%

     ▁▃▅▆████▇▆▅▄▃▂▁▁▁▂▂▂▂▂       ▁▁▁▁                         ▂
  ▅▆█████████████████████████▇█▇▇███████▇▇▇▆▅▆▆▆▆▅▆▆▅▂▄▄▅▃▄▄▆▅ █
  407 ns        Histogram: log(frequency) by time       466 ns <

 Memory estimate: 576 bytes, allocs estimate: 3.

julia> p = rand(20);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 171 evaluations.
 Range (min … max):  632.795 ns …   2.940 μs  ┊ GC (min … max): 0.00% … 77.83%
 Time  (median):     647.421 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   656.425 ns ± 106.955 ns  ┊ GC (mean ± σ):  0.79% ±  3.76%

       ▂▂▄▆▇█▇▇▆▄▃▁▁▁▁▁▃▃▃▂▁▁    ▁▁▁▂                           ▂
  ▂▄▅▆▇███████████████████████████████▇▇▇▇▅▆▆▆▆▆▅▆▅▅▄▅▄▄▂▅▄▃▅▄▄ █
  633 ns        Histogram: log(frequency) by time        713 ns <

 Memory estimate: 672 bytes, allocs estimate: 3.

julia> p = rand(24);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 45 evaluations.
 Range (min … max):  908.333 ns …  10.032 μs  ┊ GC (min … max): 0.00% … 90.52%
 Time  (median):     934.244 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   945.136 ns ± 218.800 ns  ┊ GC (mean ± σ):  0.56% ±  2.21%

        ▁▄▄▅▆▇███▇▆▄▂▁    ▁▃▃▄▃▃▃▁▁                             ▂
  ▄▃▄▄▆▇██████████████▇▅▆▇██████████████▇█▇▇█▇█▇▇█▆▅▅▄▃▄▅▅▅▄▆▅▆ █
  908 ns        Histogram: log(frequency) by time       1.02 μs <

 Memory estimate: 768 bytes, allocs estimate: 3.

So the unrolled tuple version is over 100x faster (of course this is much more memory friendly)

Would that vectorize better than a fully unrolled one? Looking at the crazy g1, it looks to me that you could write that function without a ton of repetition using Base.Cartesian.@nexprs for the declaration and the loops and stuff.

Wow, well honestly a lot of that is pretty over my head here. I have never really taken SIMD seriously. I'm trying to write the @nexprs equivalent of your g1 and I have one stumper:
function g1_gen(s)
  # declare the a and b scalar variables:
  Base.Cartesian.@nexprs 17 j->(b_j = inv(s[j+1] - s[j]))
  Base.Cartesian.@nexprs 17 j->(a_j = b_j*s[j])
  # first sequence of operations for index=1:
  CC1 = CC[1]
  Base.Cartesian.@nexprs 16 j->begin
    a_j = muladd(a_j, -CC1[j], a_{j+1})
    b_j = muladd(a_j, -CC1[j], b_{j+1})
  end
  # outer loop over "index".
  # STUMPER: Is there some way to make sure that this gets
  # expanded before the "inner" @nexprs, so that {16-index_m1} expands to the
  # raw int before the inner call to @nexprs?
  Base.Cartesian.@nexprs 15 index_m1->begin
    index = index_m1+1
    # inner loop over the variable index:
    Base.Cartesian.@nexprs {16-index_m1} j->begin
      a_j = a_{j+1}*CC[index][j]
      b_j = b_{j+1}*CC[index][j]
    end
  end
  a_1/b_1
end

But if there is some trick to make that outer @nexprs expand first, then I think this will give what you want. I'm sure it will need to be workshopped plenty, but just playing around for fun. I tried putting a @macroexpand on the outer @nexprs, but that didn't change anything.

Thanks for the idea 😄 I think this can go into an example of the what the heck / or amazing examples of julia but
function g1(s)
    @nexprs 17 i -> b_{i} = 1 / (s[i + 1] - s[i])
    @nexprs 17 i -> a_{i} = s[i] * b_{i}
    @nexprs 16 index -> (@nexprs 17-index i -> (a_{i}, b_{i}) = (muladd(a_{i}, -CC[index][i], a_{i+1}), muladd(b_{i}, -CC[index][i], b_{i+1})))
    return a_1 / b_1
end

I believe that is correct (at least it returns correct answer ) and what you were looking for in like three lines of code!
julia> @benchmark g1($(0.11274150564997776, 0.11475646021904119, 0.11370868384312821, 0.11451023277070166, 0.11369906525599734, 0.11472113632452476, 0.11317926911257481, 0.11588910073757686, 0.11045137194340607, 0.12271888810305538, 0.09198318397942494, 0.17666004884002684, -0.07776136204113551, 0.7501986007978472, -2.150973108989948, 8.739300196625988, -34.86022928456346,<PHONE_NUMBER>7542904))
BenchmarkTools.Trial: 10000 samples with 991 evaluations.
 Range (min … max):  41.835 ns … 73.957 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     41.961 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.132 ns ±  1.291 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄█▆▃▂                                                       ▁
  █████▅▁▁▁▁▁▁▁▄▁▁▁▃▃▁▃▃▄▄▄▄▃▅▄▁▁▁▄▃▁▁▄▁▅███▅▅▇▄▁▁▁▁▃▄▁▅▃▁▃▃▄ █
  41.8 ns      Histogram: log(frequency) by time        46 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

I think we could use a generated function as well based on the NTuple length of N.
@oscardssmith I'm curious what you feel like including code like that in the package considering compile times / package load times / cache pressure? The generated llvm code is immense and how does that affect other code who may call besselk amongst other calls within a larger function. It is definitely very very fast in these microbenchmarks...

in terms of speed, the generated function + nexprs approach will be really good. and I don't think the cache pressure/load times will be much worse than the vectorized evalpolys (which also unroll). The one thing you could do that would help is generating the n^2 nexprs as vectorized muladds in the first place (since doing so will reduce the number of instructions by 4x/8x).

I can't believe it. I just had to stop trying to be clever with the {} or the [] and just write the raw difference! Ugh. I still feel like the curly brace definition should have worked based on the docs for @nexprs. But anyways, your version looks pretty slick! And I'm glad that the perf aspect is appealing.

Haha ya one thing about these type of expressions is they aren't forgiving :) but ya this is really elegant. I like it a lot and it's very fast. I would think it would be good to move forward with this type of approach. I will work on Oscar's suggestion a little bit to reduce the number of instructions but here is a fully general version of Levin's sequence transformation (which works for any N). And also because we are unrolling it like this we can just put the constant factor within the function at no cost which helps it be fully general.

@inline levin_scale(B, n, k) = -(B + n) * (B + n + k)^(k - 1.0) / (B + n + k + 1.0)^k

@generated function levin_transform(s::NTuple{N, Float64}) where N
    len = N-1
    len2 = N-2
    :(
        begin
            @nexprs $len i -> b_{i} = 1 / (s[i + 1] - s[i])
            @nexprs $len i -> a_{i} = s[i] * b_{i}
            @nexprs $len2 k -> (@nexprs ($len-k) i -> (a_{i}, b_{i}) = (muladd(a_{i}, Levin_scale(1.0, i, k-1), a_{i+1}), muladd(b_{i}, Levin_scale(1.0, i, k-1), b_{i+1})))
            return a_1 / b_1
        end
    )
end

One thing to note is we aren't quite done yet is that this is a sequence transform so we need to generate efficiently the sequence (cumsum) of our series definition. Shouldn't be too time consuming I think the straightforward approach of going from series definition to sequence will work just fine (we will also have a better estimate of our weight function here doing from the series definition).

Ok here is an actual version that we can test and improve to compute besselk using the Levin's sequence transformation. It computes the partial sums and weights on the fly and can be run for any length of order.
This version works with SIMDMath.jl now to explicitly vectorize all these calculations. O and I have a version that works for complex numbers!

using Base.Cartesian
using SIMDMath: VE, muladd, LVec, FloatTypes

@inline levin_scale(B::T, n, k) where T = -(B + n) * (B + n + k)^(k - one(T)) / (B + n + k + one(T))^k

# compute the scaled modified bessel function using the Levin-u transformation
# v = 0.2; x = 4.3
# besselkx(v, x, Val(19))
@generated function besselkx_levin(v, x::T, ::Val{N}) where {T <: FloatTypes, N}
    len = N-1
    len2 = N-2
    :(
        begin
            # generate a sequence of partial sums and weights of the asymptotic expansion for large argument besselk
            l = let (out, term, fv2) = (one(typeof(x)), one(typeof(x)), 4*v^2)
                @ntuple $N i -> begin
                    offset = muladd(2, i, -1)
                    term *= muladd(offset, -offset, fv2) / (8 * x * i)
                    out += term
                    # this check is relatively cheap but if the normal sum converges then we should return value
                    abs(term) <= eps(T) && return out * sqrt(π / 2x)
                    # need to return the weights and weighted partial sums
                    invterm = inv(term)
                    LVec{2, T}((out * invterm, invterm))
                end
            end
            # this is required because you can't set the index of a tuple
            @nexprs $len i -> a_{i} = l[i]
            @nexprs $len2 k -> (@nexprs ($len-k) i -> a_{i} = muladd(a_{i}, levin_scale(one(T), i, k-1), a_{i+1}))
            return (a_1[1].value / a_1[2].value) * sqrt(π / 2x)
        end
    )
end

@generated function besselkx_levin(v, x::Complex{T}, ::Val{N}) where {T <: FloatTypes, N}
    len = N-1
    len2 = N-2
    :(
        begin
            # generate a sequence of partial sums and weights of the asymptotic expansion for large argument besselk
            l = let (out, term, fv2) = (one(typeof(x)), one(typeof(x)), 4*v^2)
                @ntuple $N i -> begin
                    offset = muladd(2, i, -1)
                    term *= @fastmath muladd(offset, -offset, fv2) / (8 * x * i)
                    out += term
                    # this check is relatively cheap but if the normal sum converges then we should return value
                    (abs(real(term)) <= eps(T) && abs(imag(term)) <= eps(T)) && return out * sqrt(π / 2x)
                    # need to return the weights and weighted partial sums
                    invterm = @fastmath inv(term)
                    LVec{4, T}((reim(out * invterm)..., reim(invterm)...))
                end
            end
            # this is required because you can't set the index of a tuple
            @nexprs $len i -> a_{i} = l[i]
            @nexprs $len2 k -> (@nexprs ($len-k) i -> a_{i} = muladd(a_{i}, levin_scale(one(T), i, k-1), a_{i+1}))
            return (complex(a_1[1].value, a_1[2].value) / complex(a_1[3].value, a_1[4].value)) * sqrt(π / 2x)
        end
    )
end

besselk_levin(v, x::T, ::Val{N}) where {T <: Union{Complex{FloatTypes}, FloatTypes}, N} = exp(-x) * besselkx_levin(v, x, Val(N))

I've added some comments in there about what it is actually doing (though still somewhat terse) but I think this is ready to be workshopped a little bit. The real version will generate <2 x double> instructions where the complex will generate <4 x double> instructions. The main challenge with the complex code is the runtime is dominated by the divisions in the sequence generation not in the transformation.
Some examples and benchmarks.
julia> v = 0.2
0.2

julia> x = 4.3
4.3

julia> besselkx_levin(v, x, Val(16))
0.5911957440242439

julia> using SpecialFunctions

julia> besselkx(v, x)
0.591195744024244

julia> using BenchmarkTools

julia> @benchmark besselkx_levin($v, $x, Val(16))
BenchmarkTools.Trial: 10000 samples with 987 evaluations.
 Range (min … max):  50.411 ns … 90.210 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     50.954 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   51.462 ns ±  4.162 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▅ ▁                                                        ▁
  ██▁██▄▄▁▄▁▁▁▁▁▁▃▁▃▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▃▁▁▁▁▃█ █
  50.4 ns      Histogram: log(frequency) by time      86.3 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And now for complex numbers!!!
julia> x = 4.3 + 1.1im
4.3 + 1.1im

julia> besselkx_levin(v, x, Val(16))
0.5783561826343409 - 0.06993526928128495im

julia> besselkx(v, x)
0.578356182634341 - 0.0699352692812849im

So you might think that the complex version will be much much slower considering all the divisions and stuff but because we are able to vectorize it pretty well it is only about 3x slower than the real version (I'm sure I can play around with the sequence generation to get a bit faster).
julia> @benchmark besselkx_levin($v, $x, Val(16))
BenchmarkTools.Trial: 10000 samples with 803 evaluations.
 Range (min … max):  155.110 ns … 180.239 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     155.340 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   155.431 ns ±   0.667 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▄▆███▆▄▁                                               ▁▁   ▂
  ▄█████████▃▁▃▃▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▄▇█████▇ █
  155 ns        Histogram: log(frequency) by time        158 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Now that I'm looking at the code more we can just combine the real and complex version and have a branch because of the type information that will be eliminated at runtime for the versions.

Okay, not sure if this is the right place for it, but I've been playing with ForwardDiff AD on this routine and I have found one issue: half-integer orders. The reason for this is because the uniform and argument-based asymptotic expansions actually terminate after a finite number of terms for v+1/2 an integer. But of course that doesn't mean partial information doesn't continue. So the issue is with
[...]
abs(term) <= eps(T) && return out * sqrt(π / 2x)
[...]

But if you just comment that out, then the inv(term) now blows up and you end up getting a NaN out once all is said an done.
My next naive thought is to just stop updating the things stored in the tuple when the series has converged:
[...]
if abs(term) > eps(T)
  out += term
  invterm = inv(term)
end
[...]

And just by putting that little branch there now the tuple values will just stay constant after term is zero.
BUT, the issue that still comes up is with the actual Levin transform stuff. And this is the part that I don't understand super well, but each of those steps starts at the first index, so even though the series is converged, you still get different numbers if you compare besselkx_levin(3.5, 4.0, Val(16)) with besselkx_levin(3.5, 4.0, Val(20)).
I don't think this is impossible, but just an update because it will probably be relevant in other places too. I think we just need to find a way to let terms be zero in peace. I'll keep tinkering with this for now though, unless your rewrite of the Levin code will fix this for free you think.

Hmmmm... I actually don't think this should be an issue with the Levin transform in itself but the generation of the partial sums. For instance if we look at the expansion in NIST 10.40.2
function _besselk_large_argument(v, x::T) where T
    MaxIter = 50
    S = eltype(x)
    v, x = S(v), S(x)

    fv2 = 4 * v^2
    term = one(S)
    res = term
    s = term
    for i in 1:MaxIter
        @show res
        offset = muladd(2, i, -1)
        term *= muladd(offset, -offset, fv2) / (8 * x * i)
        res = muladd(term, s, res)
        #abs(term) <= eps(T) && break
    end
    return res
end

These terms just generate ones. Which when we plug in the weight function this becomes zeros for the weights. Now zeros for the weights will just generate NaNs really quickly.  But what we need to have happen is to differentiate to the form given in 10.40.8 which I showed in the other issue before calculating the weights..... So to me if ForwardDiff can work with the regular asymptotic expansion it should also work here. Our partial sums we generate our just not right. This seems to be an issue with how I've set the code up and not in the actual transform itself? Will look more at this.

For sure, I'm with you about that---I just thought there was something about the Levin transform that required keeping that second invterm. I'll also poke around a bit more with this and I'm sure we can sort it out.

Could you also post code that produces wrong results? Like does it produce zeros or the wrong results and is this just for v=1/2 ? I just tried with enzyme and it seemed reasonable?
 julia> autodiff(Forward, v-> besselkx_levin(v, 12.0, Val(16)), Duplicated, Duplicated(0.5, 1.0))
(0.3618006272791338, 0.015075026136630575)

julia> besselkx(0.5, 12.0)
0.3618006272791339

julia> autodiff(Forward, v-> besselkx_levin(v, 12.0, Val(16)), Duplicated, Duplicated(1.5, 1.0))
(0.39195067955239493, 0.04710945667697055)

julia> besselkx(1.5, 12.0)
0.3919506795523951

I'll write more on this tomorrow, but unfortunately the answers that seem reasonable are often wrong for this edge case. Consider:
using FiniteDifferences, Bessels
const FDM = central_fdm(10,1)
function ref_dbesselkx_dv(v, x)
  f2 = Base.Fix2(Bessels.besselkx, x)
  FDM(f2, v)
end

ref_dbesselkx_dv(0.5, 12.0) # ~ 0.0145
ref_dbesselkx_dv(1.5, 12.0) # ~ 0.0470

So what you're seeing is not catastrophically inaccurate, but also only correct to a couple significant digits. Here is the version of your code that I am using to play with ForwardDiff (really just cosmetic modifications to fit my own formatting tastes):
@inline function levin_scale(B::T, n, k) where T
  -(B + n) * (B + n + k)^(k - one(T)) / (B + n + k + one(T))^k
end

@generated function _besselkx_levin(v::V, x::T, ::Val{N}) where {V, T, N}
  len  = N-1
  len2 = N-2
  outT = promote_type(V, T)
  quote
    begin
      l = let (out, term, fv2, eightx) = (one($outT), one($outT), 4*v^2, 8*x)
        Base.Cartesian.@ntuple $N i -> begin
          offset = muladd(2, i, -1)
          term  *= muladd(offset, -offset, fv2) / (eightx*i)
          out   += term
          abs(term) <= 1e-15 && return out * sqrt(π / 2x)
          invterm = inv(term)
          (out*invterm, invterm)
        end
      end
      Base.Cartesian.@nexprs $len i -> a_{i} = l[i]
      Base.Cartesian.@nexprs $len2 k -> begin
        Base.Cartesian.@nexprs ($len-k) i -> begin
          sc  = levin_scale(one(T), i, k-1)
          a_i = (muladd(a_{i}[1], sc, a_{i+1}[1]),
                 muladd(a_{i}[2], sc, a_{i+1}[2]))
        end
      end
      return $outT((a_1[1]/a_1[2]) * sqrt(π / 2x))
    end
  end
end

function our_dbesselkx_dv(v, x)
  f2 = _v -> _besselkx_levin(_v, x, Val(16))
  ForwardDiff.derivative(f2, v)
end

And again you see:
julia> ref_vs_ours(v, x) = (ref_dbesselkx_dv(v, x) - our_dbesselkx_dv(v, x))/abs(ref_dbesselkx_dv(v, x))
ref_vs_ours (generic function with 1 method)

julia> ref_vs_ours(1.5, 12.0)
-0.0020240634451386256

julia> ref_vs_ours(1.501, 12.0)
-2.9332885207854924e-13

The reason for this is the early return in the Levin code, and so if the series happens to terminate exactly (which occurs whenever v+1/2 is an integer), you hit that early return and you effectively differentiate just the asymptotic expansion itself, and with only ceil(v) many terms. And the reason THAT is wrong is because you could write $\mathcal{K}_{\nu}(x)$ as
$ \mathcal{K}{\nu}(x) = \sum{j=1}^{\text{ceil}(\nu)} [...] + R_{\nu}(x) $
where $R$ is some remainder term that is zero when $\nu + 1/2$ is an integer. But $\partial_{\nu} R$ is certainly not zero even when $\nu + 1/2$ is an integer, and so this early return really cuts off an entire piece of the function.
Anyways, this comment is getting long, sorry. I do think this can be addressed, but it will probably in the end involve finding a way to accumulate coefficients even when the not-derivative part is 0, if you know what I mean.

So you are saying we should just use finite differences :smile: This is an interesting problem (and at least we seem to be tackling toughest case first. I just want to take one step back to make sure I see the problem. If we look at the regular asymptotic expansion
function besselk_large_argument(v, x::T) where T
    a = exp(-x / 2)
    coef = a * sqrt(pi / 2x)
    return T(_besselk_large_argument(v, x) * coef * a)
end
function _besselk_large_argument(v, x::T) where T
    MaxIter = 50
    S = eltype(x)
    v, x = S(v), S(x)

    fv2 = 4 * v^2
    term = one(S)
    res = term
    s = term
    for i in 1:MaxIter
        offset = muladd(2, i, -1)
        term *= muladd(offset, -offset, fv2) / (8 * x * i)
        res = muladd(term, s, res)
        #abs(term) <= eps(T) && break
    end
    return res
end

And I modify the code to just run for 50 terms no matter what we can get good derivatives? At least that seems to match ok.
julia> autodiff(Forward, v->besselk_large_argument(v, 55.0), Duplicated, Duplicated(0.5, 1.0))[2]
1.97876284377238e-27

julia> dvbesselk(ArbFloat(0.5), ArbFloat(55.0))
1.97876284377238030106336346605515496107532396354231328952625258552628280879127918412731333961880277672326758481184840407518354008994240482915254186001225691837201407244641556100665654134700311440503454914579578322988897030296403705253013030240425314170276324542186929853440167150195734242472842029798051876839218023694971818866179030774694381722849809203226358554462257075547987853989971279198759455911150377220806768520043019343563369383924297628343713389132011187097644750889728058214066854106251841856730816010027578080074991956856259721740611201807109150298939419064511338703147526835965006e-27

So that is right to every digit. It does seem like the derivative is more slowly converging though so need to be cautious with that termination criteria?
So the next question then is how does that translate to the sequence transform.... I'd be interested to see if we can separate out the derivative partial sums from the actual transform itself to see if we can feed in the right partial sums if that helps converge the sum. One thing to keep in mind is that the sequence transforms only work well for diverging alternating series. It looks like the derivative will satisfy but might have slightly different convergent characteristics?

Ha, I would definitely not suggest just doing finite differences. But yeah, what you're describing about just having a bunch of terms to compensate is exactly what we do in BesselK.jl. The issue is just that for the real intermediate cases like (v,x) = (3.5, 5.0), you will just need an absolute ton of terms. I think the series transform things are super interesting there, and I don't quite understand the Miller method enough to comment on any of it. But I think this strategy is super appealing and fixable. In fact, if we could have a sequence transform that didn't involve a term in S_{n+1} - S_n in the denominator that will be zero when the series converges in finitely many steps (that are only knowable at runtime), I think we'd be good to go. We just need to think about it for a little while I think. I'm confident that there will be a way to deal with this and keep the good tricks that give us such better accuracy per iteration in the intermediate ranges.

Love this problem :smile:! Super interesting. The good news though is the sequence transform could drastically reduce the number of terms we need so I think this is still very promising especially if the derivative takes even more terms.
With that said I think this might be just an issue with how I originally coded the first problem If we take this all the way back (I'm just going to put the general code here).
function g(s)
           n = length(s) - 1
           l = similar(s)
           N = similar(l)
           D = similar(N)

           for k in 1:n
               w = s[k+1] - s[k]
               N[k], D[k] = s[k] / w, 1.0 / w
               kp1 = k + 1.0
               kdkp1 = k/kp1
               kdkp1jm2 = 1/kdkp1
               for j in 1:k-1
                   cst = -(kp1-j)*kdkp1jm2/kp1
                   N[k-j] = muladd(cst,N[k-j],N[k-j+1])
                   D[k-j] = muladd(cst,D[k-j],D[k-j+1])
                   kdkp1jm2 *= kdkp1
               end
               l[k] = N[1]/D[1]
           end
           return l
       end
function _besselk_large_argument(v, x::T, iter) where T
           MaxIter = 50
           S = eltype(x)
           v, x = S(v), S(x)

           fv2 = 4 * v^2
           term = one(S)
           res = term
           s = term
           for i in 1:iter
               offset = muladd(2, i, -1)
               term *= muladd(offset, -offset, fv2) / (8 * x * i)
               res = muladd(term, s, res)
               #abs(term) <= eps(T) && break
           end
           return res  * sqrt(pi / 2x)
       end

Now what we need to realize is that the sequence transform will be able to handle any sequence but we don't necessarily need to take the derivative of the sequence transform we only need to take the derivative of the sequence then feed that sequence into the transform. I'm hoping that makes sense.
So essentially this is all we need to do...

v = 0.5; x = 5.0
# generate our partial sums but using the derivative
out = [autodiff(Forward, v->_besselk_large_argument(v, x, iter), Duplicated, Duplicated(v, 1.0))[2] for iter in 1:30];

# do the sequence transform on those derivative partial sums
julia> g(out)
30-element Vector{Float64}:
 0.05604991216397929
 0.051379086150314356
 0.05131596687985941
 0.0513217030639802
 0.051321108923874925
 0.051321175118541326
 0.05132116743538506
 0.051321168343865145
 0.051321168236451525
 0.05132116824892151
 0.051321168247529726
 0.051321168247674714
 0.051321168247661336
 0.051321168247662294
 0.05132116824766225
 0.05132116824766225
 0.05132116824766225
 0.05132116824766225
 0.051321168247662266
 0.05132116824766226
 0.05132116824766225
 0.05132116824766226
 0.05132116824766225
 0.05132116824766226
 0.05132116824766225
 0.051321168247662266
 0.05132116824766226
 0.05132116824766225
 0.051321168247662266

julia> dvbesselk(ArbFloat(v), ArbFloat(x))[2] * exp(ArbFloat(x))
0.05132116824766225501961245911721947725432532303036730970326023167364314434751562772143246443214094748722

So as you can see this only needs 14 terms to give us a derivative accurate to every digit!

Wow, okay, very nice! So we really just need to somehow get those derivatives evaluated before the acceleration computes inv(w). That is super doable.

Yes - I think there's two things about this problem.

Separate the code for the partial sums from the sequence transformation. Use autodiff on the partial sums then feed those into the sequence transform. The Levin function can then just be a general function for any sequence transformation that accepts partical sums sn and corresponding weights wn.

How do we have a stopping criteria for: a. how many terms to use in the sequence of partial sums or b. how many terms to use in the iterative normal asymptotic expansion that converges for both the scalar evaluation and derivative.

To me the second point has to be some general problem with AD itself for iterative processes. I'm wondering if this is already solved? I'm not as familiar with AD honestly but the problem to me could be generally stated as for an iterative sum (with i terms) to what extent fi(x) converges to f(x) which we detect with our stopping criteria during the recursive sum implies that fi'(x) converges to f'(x). Ideally, AD would be able to check for both conditions.... but it seems like it just terminates when the first condition is true which generates awful estimates of the second condition.
This issue has to arise frequently? As I would mostly expect the sum of derivatives to be slower converging (at least require a couple more higher order terms). Is it possible to have two separate stopping criteria when AD is used?

Can the Implicit Function theorem be used here?

Okay, so after looking at the Enzyme docs for a while, the solution of compartmentalizing the Levin transformation from the series accumulation actually has a reasonably elegant solution. If I go back to an earlier version of the code that doesn't also collect invterm, I can do the Levin transform directly on the derivative of the series values like so:
using Bessels, Enzyme
using Enzyme.EnzymeRules
using Base.Cartesian: @nexprs, @ntuple

@inline levin_scale(B::T, n, k) where T = -(B + n + k) * (B + n + k - 1) / ((B + n + 2k) * (B + n + 2k - 1))

@generated function levin_transform(s::NTuple{N, Float64}) where N
    len = N-1
    len2 = N-2
    :(
        begin
            @nexprs $len i -> b_{i} = 1 / (s[i + 1] - s[i])
            @nexprs $len i -> a_{i} = s[i] * b_{i}
            @nexprs $len2 k -> (@nexprs ($len-k) i -> (a_{i}, b_{i}) = (muladd(a_{i}, levin_scale(1.0, i, k-1), a_{i+1}), muladd(b_{i}, levin_scale(1.0, i, k-1), b_{i+1})))
            return a_1 / b_1
        end
    )
end

# Here is the sauce: a manual rule that separates does the levin transform for
# the values and dvalues separately. This is because if v+1/2 is an integer, the
# sequences of _values_ will be exactly converged and you'll have a divide by
# zero NaN problem, but the _dvalues_ will NOT do that, and you want all those
# terms for the sake of accuracy.
function EnzymeRules.forward(func::Const{typeof(levin_transform)},
                             ::Type{<:Duplicated},
                             s::Duplicated)
  sv  = s.val
  dv  = s.dval
  N   = length(sv)
  ls  = (sv[N-1] == sv[N]) ? sv[N] : levin_transform(sv)
  dls = (dv[N-1] == dv[N]) ? dv[N] : levin_transform(s.dval)
  Duplicated(ls, dls)
end

@generated function besselkx_as(v, x::T, ::Val{N}) where {T, N}
  quote
    # generate a sequence of partial sums and weights of the asymptotic
    # expansion for large argument besselk
    tup = let (out, term, fv2) = (one(typeof(x)), one(typeof(x)), 4*v^2)
      @ntuple $N i -> begin
        offset  = muladd(2, i, -1)
        term   *= muladd(offset, -offset, fv2) / (8 * x * i)
        out    += term
        out
      end
    end
  end
end

@generated function besselkx_levin(v, x::T, ::Val{N}) where {T, N}
  quote
    tmp  = besselkx_as(v, x, Val($N))
    tmp2 = levin_transform(tmp)
    tmp2*sqrt(π / 2x)
  end
end

# This allocates, so there may still be something to chase down.
dbeskx_dv(v, x) = autodiff(Forward, _v->besselkx_levin(_v, x, Val(20)),
                           Duplicated, Duplicated(v, 1.0))[2]

#=
# A reference so that you can check that the half-integer derivatives work.
#
using FiniteDifferences
dbeskx_dv_ref(v, x) = central_fdm(10,1)(_v->besselkx(_v, x), v)
=#

And you can check that the derivatives agree with the finite diff incredibly well even for v=1.5, v=2.5, and so on. This seems like a complete win to me.
So, @heltonmc: at some point, can you elaborate to me on why you have the invterm stuff in the current besselkx_levin? We do get the right answers without it, so does it somehow significantly speed up the actual Levin transformation code?

This is awesome. Should we maybe split the particular Enzyme stuff into a separate issue to keep track of?
Just to keep everyone in the loop I think we are going to have to go the route of custom rules for some functions. That might be only way to get great performance and accuracy. The plan was to perhaps have this as a separate module that lives within Bessels and use a weak dep. and package extension or perhaps move it into it's own separate package like BesselsAD that focuses on being compatible with perhaps a larger AD framework.
@cgeoga I just pushed recent changes to the besselk_levin code over at #95. I'll paste the new functions below.

@generated function besselkx_levin(v, x::T, ::Val{N}) where {T <: FloatTypes, N}
    :(
        begin
            s = zero(T)
            t = one(T)
            l = @ntuple $N i -> begin
                    s += t
                    t *= (4*v^2 - (2i - 1)^2) / (8 * x * i)
                    invterm = inv(t)
                    Vec{2, T}((s * invterm, invterm))
                end
            return levin_transform(l) * sqrt(π / 2x)
        end
    )
end

@generated function besselkx_levin(v, x::Complex{T}, ::Val{N}) where {T <: FloatTypes, N}
    :(
        begin
            s = zero(T)
            t = one(typeof(x))
            t2 = t
            a = @fastmath inv(8*x)
            a2 = 8*x

            l =  @ntuple $N i -> begin
                    s += t
                    b = (4*v^2 - (2i - 1)^2) / i
                    t *= a * b
                    t2 *= a2 / b
                    Vec{4, T}((reim(s * t2)..., reim(t2)...))
                end
            return levin_transform(l) * sqrt(π / 2x)
        end
    )
end

These functions take advantage of constant propagation. They were a little faster but was just trying to improve the code clarity a little bit. I also refactored the complex code. The issue with the complex function is you need to desperately avoid any complex divisions so the different accumulators are to make sure we only perform real divisions or at worst a complex divided by a real and always avoid a real divided by complex or complex divided by complex. Even doing like ten extra multiplies is worth this cost. If you could try it with this new code that would be great on that other branch!

Very nice. Now that I understand how to do the levin_transform trickery I think I can write rules that work with those output Vec{N,T} objects and stuff. Ideally, I can just write a few of those manual forward rules for EnzymeRules.forward(::Const{typeof(levin_transform)}, [...] and we won't have to touch anything else. I don't expect any issues that can't be tinkered away.
Also, good point about a separate issue. I'll start a new one and use that method as the top post on it when it is done.


Create a Backup and Recovery process in an AI Farm Environment

An AI farm environment is a complex infrastructure that requires careful planning for backup and recovery. In this environment, it is important to have a reliable backup and recovery strategy to protect critical data and ensure business continuity.
One recommended strategy is to make a backup on the local host directory of each Microserver, and then make a remote backup of the microserver backup on the GPU Server. This approach provides an added layer of protection against data loss or corruption.
Additionally, it is important to regularly test the backup and recovery process to ensure that it is working effectively. By following these best practices, organizations can ensure that their AI farm environment is well-protected and prepared for any potential disasters.

Completed

Completed


Use contenteditable for editing text

First step to address #1126.
I tested it in Chrome and FF on MacOS, would be very helpful to get folks to test other browser/OS combos.
For the future, I'm planning to copy as much as I can from Quip. That means creating new elements inside the contenteditable after user selects a range and chooses a color.
Onwards and forwards.

I'm sure you're aware how hard it is to make contenteditable behave identically on all browsers (or to at least account for it). There's been multiple posts written about it (just to pick the ones I have bookmarked — [1], [2], [3]). I've had to deal with the problems myself when we were still using it, which is why we've ultimately moved to textarea.
So in short, making contenteditable work will take a great deal of work and you'll spend more time fixing incompatibilities reinventing the wheel than doing any actual development.
I think it'd be easier to pick some existing light-weight editor and build on top of that instead.

I think those problems with contenteditable are more likely to show up once we get to the editing, but it does look like a mine field.
Have you looked at which editor you'd use? Quill comes up first on Google and looks decent: https://quilljs.com/
I know at FB we devolped Draft https://draftjs.org/ but it is now considered somewhat badly designed/implemented and we're rewriting it internally.
cc @vjeux if you have any ideas.

This is probably quite outdated by now but may still be valuable in this conversation. :)

Quill sits at 45kB, although maybe our build could shave some of it off - that seems ok given Excalidraw is somewhere at 1MB+.

@kbariotis Thanks, that's awesome! Looks like good learnings from the approach in that PR. I will look into this when I get some free time.