ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Add Integrate.OnOpenInterval

uap-experimental
diluculo 7 years ago
parent
83ab6ebd97
commit
92ae8e358a
3 changed files with 153 additions and 34 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 77
      .paket/Paket.Restore.targets
  2. 38
      src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
  3. 72
      src/Numerics/Integrate.cs

77
.paket/Paket.Restore.targets
View File

@ -5,6 +5,11 @@
<!-- make MSBuild track this file for incremental builds. -->
<!-- ref https://blogs.msdn.microsoft.com/msbuild/2005/09/26/how-to-ensure-changes-to-a-custom-target-file-prompt-a-rebuild/ -->
<MSBuildAllProjects>$(MSBuildAllProjects);$(MSBuildThisFileFullPath)</MSBuildAllProjects>
<DetectedMSBuildVersion>$(MSBuildVersion)</DetectedMSBuildVersion>
<DetectedMSBuildVersion Condition="'$(MSBuildVersion)' == ''">15.0.0</DetectedMSBuildVersion>
<MSBuildSupportsHashing>false</MSBuildSupportsHashing>
<MSBuildSupportsHashing Condition=" '$(DetectedMSBuildVersion)' &gt; '15.8.0' ">true</MSBuildSupportsHashing>
<!-- Mark that this target file has been loaded. -->
<IsPaketRestoreTargetsFileLoaded>true</IsPaketRestoreTargetsFileLoaded>
<PaketToolsPath>$(MSBuildThisFileDirectory)</PaketToolsPath>
@ -59,7 +64,7 @@
<!-- Disable automagic references for F# dotnet sdk -->
<!-- This will not do anything for other project types -->
<!-- see https://github.com/fsharp/fslang-design/blob/master/RFCs/FS-1032-fsharp-in-dotnet-sdk.md -->
<!-- see https://github.com/fsharp/fslang-design/blob/master/tooling/FST-1002-fsharp-in-dotnet-sdk.md -->
<DisableImplicitFSharpCoreReference>true</DisableImplicitFSharpCoreReference>
<DisableImplicitSystemValueTupleReference>true</DisableImplicitSystemValueTupleReference>
@ -73,37 +78,47 @@
<MSBuild Projects="$(PaketToolsPath)paket.bootstrapper.proj" Targets="Restore" />
</Target>
<!-- Official workaround for https://docs.microsoft.com/en-us/visualstudio/msbuild/getfilehash-task?view=vs-2019 -->
<UsingTask TaskName="Microsoft.Build.Tasks.GetFileHash" AssemblyName="Microsoft.Build.Tasks.Core, Version=15.1.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(DetectedMSBuildVersion)' &lt; '16.0.360' " />
<UsingTask TaskName="Microsoft.Build.Tasks.VerifyFileHash" AssemblyName="Microsoft.Build.Tasks.Core, Version=15.1.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(DetectedMSBuildVersion)' &lt; '16.0.360' " />
<Target Name="PaketRestore" Condition="'$(PaketRestoreDisabled)' != 'True'" BeforeTargets="_GenerateDotnetCliToolReferenceSpecs;_GenerateProjectRestoreGraphPerFramework;_GenerateRestoreGraphWalkPerFramework;CollectPackageReferences" DependsOnTargets="PaketBootstrapping">
<!-- Step 1 Check if lockfile is properly restored -->
<!-- Step 1 Check if lockfile is properly restored (if the hash of the lockfile and the cache-file match) -->
<PropertyGroup>
<PaketRestoreRequired>true</PaketRestoreRequired>
<NoWarn>$(NoWarn);NU1603;NU1604;NU1605;NU1608</NoWarn>
<CacheFilesExist>false</CacheFilesExist>
<CacheFilesExist Condition=" Exists('$(PaketRestoreCacheFile)') And Exists('$(PaketLockFilePath)') ">true</CacheFilesExist>
</PropertyGroup>
<!-- Because ReadAllText is slow on osx/linux, try to find shasum and awk -->
<PropertyGroup>
<PaketRestoreCachedHasher Condition="'$(OS)' != 'Windows_NT' And '$(PaketRestoreCachedHasher)' == '' And Exists('/usr/bin/shasum') And Exists('/usr/bin/awk')">/usr/bin/shasum "$(PaketRestoreCacheFile)" | /usr/bin/awk '{ print $1 }'</PaketRestoreCachedHasher>
<PaketRestoreLockFileHasher Condition="'$(OS)' != 'Windows_NT' And '$(PaketRestoreLockFileHash)' == '' And Exists('/usr/bin/shasum') And Exists('/usr/bin/awk')">/usr/bin/shasum "$(PaketLockFilePath)" | /usr/bin/awk '{ print $1 }'</PaketRestoreLockFileHasher>
<!-- Read the hash of the lockfile -->
<GetFileHash Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(CacheFilesExist)' == 'true' " Files="$(PaketLockFilePath)" Algorithm="SHA256" HashEncoding="hex" >
<Output TaskParameter="Hash" PropertyName="PaketRestoreLockFileHash" />
</GetFileHash>
<!-- Read the hash of the cache, which is json, but a very simple key value object -->
<PropertyGroup Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(CacheFilesExist)' == 'true' ">
<PaketRestoreCachedContents>$([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)'))</PaketRestoreCachedContents>
</PropertyGroup>
<ItemGroup Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(CacheFilesExist)' == 'true' ">
<!-- Parse our simple 'paket.restore.cached' json ...-->
<PaketRestoreCachedSplitObject Include="$([System.Text.RegularExpressions.Regex]::Split(`$(PaketRestoreCachedContents)`, `{|}|,`))"></PaketRestoreCachedSplitObject>
<!-- Keep Key, Value ItemGroup-->
<PaketRestoreCachedKeyValue Include="@(PaketRestoreCachedSplitObject)"
Condition=" $([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `&quot;: &quot;`).Length) &gt; 1 ">
<Key>$([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[0].Replace(`"`, ``).Replace(` `, ``))</Key>
<Value>$([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[1].Replace(`"`, ``).Replace(` `, ``))</Value>
</PaketRestoreCachedKeyValue>
</ItemGroup>
<PropertyGroup Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(CacheFilesExist)' == 'true' ">
<!-- Retrieve the hashes we are interested in -->
<PackagesDownloadedHash Condition=" '%(PaketRestoreCachedKeyValue.Key)' == 'packagesDownloadedHash' ">%(PaketRestoreCachedKeyValue.Value)</PackagesDownloadedHash>
<ProjectsRestoredHash Condition=" '%(PaketRestoreCachedKeyValue.Key)' == 'projectsRestoredHash' ">%(PaketRestoreCachedKeyValue.Value)</ProjectsRestoredHash>
</PropertyGroup>
<!-- If shasum and awk exist get the hashes -->
<Exec StandardOutputImportance="Low" Condition=" '$(PaketRestoreCachedHasher)' != '' " Command="$(PaketRestoreCachedHasher)" ConsoleToMSBuild='true'>
<Output TaskParameter="ConsoleOutput" PropertyName="PaketRestoreCachedHash" />
</Exec>
<Exec StandardOutputImportance="Low" Condition=" '$(PaketRestoreLockFileHasher)' != '' " Command="$(PaketRestoreLockFileHasher)" ConsoleToMSBuild='true'>
<Output TaskParameter="ConsoleOutput" PropertyName="PaketRestoreLockFileHash" />
</Exec>
<!-- Debug whats going on -->
<Message Importance="low" Text="calling paket restore with targetframework=$(TargetFramework) targetframeworks=$(TargetFrameworks)" />
<PropertyGroup Condition="Exists('$(PaketRestoreCacheFile)') ">
<!-- if no hash has been done yet fall back to just reading in the files and comparing them -->
<PaketRestoreCachedHash Condition=" '$(PaketRestoreCachedHash)' == '' ">$([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)'))</PaketRestoreCachedHash>
<PaketRestoreLockFileHash Condition=" '$(PaketRestoreLockFileHash)' == '' ">$([System.IO.File]::ReadAllText('$(PaketLockFilePath)'))</PaketRestoreLockFileHash>
<PropertyGroup Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(CacheFilesExist)' == 'true' ">
<!-- If the restore file doesn't exist we need to restore, otherwise only if hashes don't match -->
<PaketRestoreRequired>true</PaketRestoreRequired>
<PaketRestoreRequired Condition=" '$(PaketRestoreLockFileHash)' == '$(PaketRestoreCachedHash)' ">false</PaketRestoreRequired>
<PaketRestoreRequired Condition=" '$(PaketRestoreLockFileHash)' == '$(ProjectsRestoredHash)' ">false</PaketRestoreRequired>
<PaketRestoreRequired Condition=" '$(PaketRestoreLockFileHash)' == '' ">true</PaketRestoreRequired>
</PropertyGroup>
@ -116,7 +131,10 @@
</PropertyGroup>
<!-- Do a global restore if required -->
<Warning Text="This version of MSBuild (we assume '$(DetectedMSBuildVersion)' or older) doesn't support GetFileHash, so paket fast restore is disabled." Condition=" '$(MSBuildSupportsHashing)' != 'true' " />
<Error Text="Stop build because of PAKET_ERROR_ON_MSBUILD_EXEC and we always call the bootstrapper" Condition=" '$(PAKET_ERROR_ON_MSBUILD_EXEC)' == 'true' AND '$(PaketBootstrapperStyle)' == 'classic' AND Exists('$(PaketBootStrapperExePath)') AND !(Exists('$(PaketExePath)'))" />
<Exec Command='$(PaketBootStrapperCommand)' Condition=" '$(PaketBootstrapperStyle)' == 'classic' AND Exists('$(PaketBootStrapperExePath)') AND !(Exists('$(PaketExePath)'))" ContinueOnError="false" />
<Error Text="Stop build because of PAKET_ERROR_ON_MSBUILD_EXEC and we need a full restore (hashes don't match)" Condition=" '$(PAKET_ERROR_ON_MSBUILD_EXEC)' == 'true' AND '$(PaketRestoreRequired)' == 'true' AND '$(PaketDisableGlobalRestore)' != 'true'" />
<Exec Command='$(PaketCommand) restore' Condition=" '$(PaketRestoreRequired)' == 'true' AND '$(PaketDisableGlobalRestore)' != 'true' " ContinueOnError="false" />
<!-- Step 2 Detect project specific changes -->
@ -126,7 +144,7 @@
<MyTargetFrameworks Condition="'$(TargetFrameworks)' != '' AND '$(TargetFramework)' == '' " Include="$(TargetFrameworks)"></MyTargetFrameworks>
<PaketResolvedFilePaths Include="@(MyTargetFrameworks -> '$(PaketIntermediateOutputPath)\$(MSBuildProjectFile).%(Identity).paket.resolved')"></PaketResolvedFilePaths>
</ItemGroup>
<Message Importance="low" Text="MyTargetFrameworks=@(MyTargetFrameworks) PaketResolvedFilePaths=@(PaketResolvedFilePaths)" />
<PropertyGroup>
<PaketReferencesCachedFilePath>$(PaketIntermediateOutputPath)\$(MSBuildProjectFile).paket.references.cached</PaketReferencesCachedFilePath>
<!-- MyProject.fsproj.paket.references has the highest precedence -->
@ -163,6 +181,7 @@
<!-- Step 3 Restore project specific stuff if required -->
<Message Condition=" '$(PaketRestoreRequired)' == 'true' " Importance="low" Text="Detected a change ('$(PaketRestoreRequiredReason)') in the project file '$(MSBuildProjectFullPath)', calling paket restore" />
<Error Text="Stop build because of PAKET_ERROR_ON_MSBUILD_EXEC and we detected a change ('$(PaketRestoreRequiredReason)') in the project file '$(MSBuildProjectFullPath)'" Condition=" '$(PAKET_ERROR_ON_MSBUILD_EXEC)' == 'true' AND '$(PaketRestoreRequired)' == 'true' " />
<Exec Command='$(PaketCommand) restore --project "$(MSBuildProjectFullPath)" --output-path "$(PaketIntermediateOutputPath)" --target-framework "$(TargetFrameworks)"' Condition=" '$(PaketRestoreRequired)' == 'true' AND '$(TargetFramework)' == '' " ContinueOnError="false" />
<Exec Command='$(PaketCommand) restore --project "$(MSBuildProjectFullPath)" --output-path "$(PaketIntermediateOutputPath)" --target-framework "$(TargetFramework)"' Condition=" '$(PaketRestoreRequired)' == 'true' AND '$(TargetFramework)' != '' " ContinueOnError="false" />
@ -224,8 +243,6 @@
<Target Name="PaketDisableDirectPack" AfterTargets="_IntermediatePack" BeforeTargets="GenerateNuspec" Condition="('$(IsPackable)' == '' Or '$(IsPackable)' == 'true') And Exists('$(PaketIntermediateOutputPath)/$(MSBuildProjectFile).references')" >
<PropertyGroup>
<ContinuePackingAfterGeneratingNuspec>false</ContinuePackingAfterGeneratingNuspec>
<DetectedMSBuildVersion>$(MSBuildVersion)</DetectedMSBuildVersion>
<DetectedMSBuildVersion Condition="$(MSBuildVersion) == ''">15.8.0</DetectedMSBuildVersion>
</PropertyGroup>
</Target>
@ -252,10 +269,12 @@
</PropertyGroup>
<ItemGroup>
<_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion).nuspec"/>
<_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion.Split(`+`)[0]).nuspec"/>
</ItemGroup>
<Exec Command='$(PaketCommand) fix-nuspecs files "@(_NuspecFiles)" project-file "$(PaketProjectFile)" ' Condition="@(_NuspecFiles) != ''" />
<Error Text="Error Because of PAKET_ERROR_ON_MSBUILD_EXEC (not calling fix-nuspecs)" Condition=" '$(PAKET_ERROR_ON_MSBUILD_EXEC)' == 'true' " />
<Exec Condition="@(_NuspecFiles) != ''" Command='$(PaketCommand) fix-nuspecs files "@(_NuspecFiles)" project-file "$(PaketProjectFile)" ' />
<Error Condition="@(_NuspecFiles) == ''" Text='Could not find nuspec files in "$(AdjustedNuspecOutputPath)" (Version: "$(PackageVersion)"), therefore we cannot call "paket fix-nuspecs" and have to error out!' />
<ConvertToAbsolutePath Condition="@(_NuspecFiles) != ''" Paths="@(_NuspecFiles)">
<Output TaskParameter="AbsolutePaths" PropertyName="NuspecFileAbsolutePath" />
@ -281,7 +300,7 @@
DevelopmentDependency="$(DevelopmentDependency)"
BuildOutputInPackage="@(_BuildOutputInPackage)"
TargetPathsToSymbols="@(_TargetPathsToSymbols)"
SymbolPackageFormat="symbols.nupkg"
SymbolPackageFormat="$(SymbolPackageFormat)"
TargetFrameworks="@(_TargetFrameworks)"
AssemblyName="$(AssemblyName)"
PackageOutputPath="$(PackageOutputAbsolutePath)"
@ -328,7 +347,7 @@
DevelopmentDependency="$(DevelopmentDependency)"
BuildOutputInPackage="@(_BuildOutputInPackage)"
TargetPathsToSymbols="@(_TargetPathsToSymbols)"
SymbolPackageFormat="symbols.nupkg"
SymbolPackageFormat="$(SymbolPackageFormat)"
TargetFrameworks="@(_TargetFrameworks)"
AssemblyName="$(AssemblyName)"
PackageOutputPath="$(PackageOutputAbsolutePath)"

38
src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
View File

@ -59,7 +59,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
{
return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y));
}
/// <summary>
/// Test Function Start point.
/// </summary>
@ -311,5 +311,41 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
Assert.AreEqual(gaussLegendre.IntervalEnd, StopA);
}
// integral_(-oo)^(oo) exp(-x^2/2) dx = sqrt(2 ¥ð)
[TestCase(double.NegativeInfinity, double.PositiveInfinity, Constants.Sqrt2Pi)]
// integral_(-oo)^(0) exp(-x^2/2) dx = sqrt(¥ð/2)
[TestCase(double.NegativeInfinity, 0, Constants.SqrtPiOver2)]
// integral_(0)^(oo exp(-x^2/2) dx = sqrt(¥ð/2)
[TestCase(0, double.PositiveInfinity, Constants.SqrtPiOver2)]
// integral_(-1)^(1) exp(-x^2/2) dx = sqrt(2 ¥ð) erf(1/sqrt(2))
[TestCase(-1, 1, 1.7112487837842976063)]
// integral_(1)^(0) exp(-x^2/2) dx = -sqrt(¥ð/2) erf(1/sqrt(2))
[TestCase(1, 0, -0.85562439189214880317)]
public void TestGaussianIntegralBySubstitution(double a, double b, double expected)
{
Assert.AreEqual(
expected,
Integrate.OnOpenInterval((x) => Math.Exp(-x * x / 2), a, b),
1e-10,
"Integral e^(-x^2 /2) from {0} to {1}", a, b);
}
// integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx
// = 1 / pi integral_(oo)^(oo) 1 / (1 + t^2) dt
// = 1
// or = 2 / pi integral_(0)^(oo) 1 / (1 + t^2) dt
// or = 2 / pi integral_(-oo)^(0) 1 / (1 + t^2) dt
[TestCase(double.NegativeInfinity, double.PositiveInfinity, 1, Constants.InvPi)]
[TestCase(0, double.PositiveInfinity, 1, Constants.TwoInvPi)]
[TestCase(double.NegativeInfinity, 0, 1, Constants.TwoInvPi)]
public void TestSincIntegralBySubstitution(double a, double b, double expected, double factor)
{
Assert.AreEqual(
expected,
factor * Integrate.OnOpenInterval((x) => 1 / (1 + x * x), a, b),
1e-10,
"Integral sin(pi*x)/(pi*x) from -oo to oo");
}
}
}

72
src/Numerics/Integrate.cs
View File

@ -45,21 +45,85 @@ namespace MathNet.Numerics
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="targetAbsoluteError">The expected relative accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double OnClosedInterval(Func<double, double> f, double intervalBegin, double intervalEnd, double targetAbsoluteError)
public static double OnClosedInterval(Func<double, double> f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError);
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth function on a closed interval.
/// Approximation of the definite integral of an analytic smooth function by substitution. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
/// <param name="f">The analytic smooth function to integrate.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="targetAbsoluteError">The expected relative accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double OnClosedInterval(Func<double, double> f, double intervalBegin, double intervalEnd)
public static double OnOpenInterval(Func<double, double> f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 1e-8);
// Reference:
// Formula used for variable subsitution from
// 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140.
// 2. quadgk.m, GNU Octave
if (intervalBegin > intervalEnd)
{
return -OnOpenInterval(f, intervalEnd, intervalBegin, targetAbsoluteError);
}
// (-oo, oo) => [-1, 1]
//
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
{
Func<double, double> u = (t) =>
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
}
// [a, oo) => [0, 1]
//
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt
// = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd))
{
Func<double, double> u = (s) =>
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
return OnClosedInterval(u, 0d, 1d, targetAbsoluteError);
}
// (-oo, b] => [-1, 0]
//
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt
// = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin))
{
Func<double, double> u = (s) =>
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
return OnClosedInterval(u, -1d, 0d, targetAbsoluteError);
}
// [a, b] => [-1, 1]
//
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else
{
Func<double, double> u = (t) =>
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
}
}
/// <summary>

Write Preview
Loading…
Cancel
Save