我正在尝试为D编程语言调试我的并行库.最近提交了一个错误报告,指出使用任务执行的某些浮点运算的低位比特在运行期间是不确定的.(如果您阅读报告,请注意并行缩减通过以确定的方式创建任务而在幕后工作.)

这似乎不是舍入模式问题,因为我尝试手动设置舍入模式.我也很确定这不是一个并发错误.该库是经过严格测试(包括传递金克斯压力测试),这个问题始终限制在低位,这甚至发生在单核的机器,其中低级别的内存模型问题是一个问题较少.浮点结果可能因调度操作的线程而异的原因还有哪些？

编辑:我在这里做一些printf调试,看起来各个任务的结果有时会在运行中有所不同.

编辑#2:以下代码以更简单的方式再现此问题.它总结了主线程中数组的术语,然后启动一个新线程来执行完全相同的函数.问题绝对不是我的库中的错误,因为这段代码甚至不使用我的库.

import std.algorithm, core.thread, std.stdio, core.stdc.fenv;

real sumRange(const(real)[] range) {
    writeln("Rounding mode:  ", fegetround);  // 0 from both threads.
    return reduce!"a + b"(range);
}

void main() {
    immutable n = 1_000_000;
    immutable delta = 1.0 / n;

    auto terms = new real[1_000_000];
    foreach(i, ref term; terms) {
        immutable x = ( i - 0.5 ) * delta;
        term = delta / ( 1.0 + x * x ) * 1;
    } …

我对Mathematica的全局优化能力有疑问.我偶然发现了与NAG工具箱相关的文本(白皮书).

现在我试图从论文中解决测试用例.正如预期的那样,Mathematica解决它的速度非常快.

n=2;
fun[x_,y_]:=10 n+(x-2)^2-10Cos[2 Pi(x-2)]+(y-2)^2-10 Cos[2 Pi(y-2)];
NMinimize[{fun[x,y],-5<= x<= 5&&-5<= y<= 5},{x,y},Method->{"RandomSearch","SearchPoints"->13}]//AbsoluteTiming

输出是

{0.0470026,{0.,{x->2.,y->2.}}}

可以看到优化例程访问的点.

{sol, pts}=Reap[NMinimize[{fun[x,y],-5<= x<= 5&&-5<= y<= 5},{x,y},Method->`{"RandomSearch","SearchPoints"->13},EvaluationMonitor:>Sow[{x,y}]]];Show[ContourPlot[fun[x,y],{x,-5.5,5.5},{y,-5.5,5.5},ColorFunction->"TemperatureMap",Contours->Function[{min,max},Range[min,max,5]],ContourLines->True,PlotRange-> All],ListPlot[pts,Frame-> True,Axes-> False,PlotRange-> All,PlotStyle-> Directive[Red,Opacity[.5],PointSize[Large]]],Graphics[Map[{Black,Opacity[.7],Arrowheads[.026],Arrow[#]}&,Partition[pts//First,2,1]],PlotRange-> {{-5.5,5.5},{-5.5,5.5}}]]`

现在我想到在更高维度上解决同样的问题.对于五个变量的问题,即使允许大量搜索点,mathematica也开始陷入局部最小陷阱.

n=5;funList[x_?ListQ]:=Block[{i,symval,rule},
i=Table[ToExpression["x$"<>ToString[j]],{j,1,n}];symval=10  n+Sum[(i[[k]]-2)^2-10Cos[2Pi(i[[k]]-2)],{k,1,n}];rule=MapThread[(#1-> #2)&,{i,x}];symval/.rule]val=Table[RandomReal[{-5,5}],{i,1,n}];vars=Table[ToExpression["x$"<>ToString[j]],{j,1,n}];cons=Table[-5<=ToExpression["x$"<>ToString[j]]<= 5,{j,1,n}]/.List-> And;NMinimize[{funList[vars],cons},vars,Method->{"RandomSearch","SearchPoints"->4013}]//AbsoluteTiming

输出不是我们想要看到的.在我的core2duo机器上花了49秒,仍然是当地的最低限度.

{48.5157750,{1.98992,{x$1->2.,x$2->2.,x$3->2.,x$4->2.99496,x$5->1.00504}}}

然后尝试了SimulatedAnealing 100000次迭代.

NMinimize[{funList[vars],cons},vars,Method->"SimulatedAnnealing",MaxIterations->100000]//AbsoluteTiming

产量仍然不合适.

{111.0733530,{0.994959,{x$1->2.,x$2->2.99496,x$3->2.,x$4->2.,x$5->2.}}}

现在,Mathematica有一个名为Minimize的精确优化算法.正如预期的那样,必须在实用性方面失败,但随着问题规模的增加,它会很快失败.

n=3;funList[x_?ListQ]:=Block[{i,symval,rule},i=Table[ToExpression["x$"<>ToString[j]],{j,1,n}];symval=10 n+Sum[(i[[k]]-2)^2-10Cos[2 Pi(i[[k]]-2)],{k,1,n}];rule=MapThread[(#1-> #2)&,{i,x}];symval/.rule]val=Table[RandomReal[{-5,5}],{i,1,n}];vars=Table[ToExpression["x$"<>ToString[j]],{j,1,n}];cons=Table[-5<=ToExpression["x$"<>ToString[j]]<= 5,{j,1,n}]/.List-> And;Minimize[{funList[vars],cons},vars]//AbsoluteTiming

输出完全没问题.

{5.3593065,{0,{x$1->2,x$2->2,x$3->2}}}

但是如果用n = 4进一步改变问题大小,你会看到结果.解决方案在我的笔记本中不会出现很长时间.

现在问题很简单,这里的任何人都认为有一种方法可以在Mathematica中有效地数值解决这个问题以获得更高维度的案例吗？让我们分享我们的想法和经验.但是应该记住,它是一个基准的非线性全局优化问题.大多数数字根寻找/最小化算法通常搜索局部最小值.

BR

P

我试图用伪谱方案数值求解Swift-Hohenberg方程http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation,其中线性项在傅立叶空间中隐式处理,而在现实空间中评估非线性.简单的Euler方案用于时间积分.
我的问题是我提出的Matlab代码完美无缺,而依赖FFTW进行傅里叶变换的C++代码变得不稳定,经过几千个步骤后就会发散.我已经跟踪了非线性项的处理方式(参见C++代码中的注释).如果我只使用Phi的真实部分,就会发生不稳定.然而,由于数值舍入误差,Phi应该只有一个可忽略的虚部,并且Matlab正在做类似的事情,保持Phi纯粹真实.在Octave下,Matlab代码也运行良好.初始条件可能类似于
R=0.02*(rand(256,256)-0.5);
Matlab(小振幅波动).

TLDR;

为什么这些代码片段做了不同的事情？具体来说,我如何使C++代码的工作方式与Matlab版本相同？

编辑1:

为了完整起见,我使用FFTW提供的R2C/C2R功能添加了代码.有关详细信息,请参阅http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html(我希望我的数据布局正确).此代码始终显示约3100个时间步后的不稳定性.如果我将dt减少到例如0.01,则会发生10次.

使用复杂DFT的C++代码

#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>

int main() {

const int N=256, nSteps=10000;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*N*sizeof(double));

// complex arrays
fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *NPhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));

// plans for Fourier transforms
fftw_plan phiPlan=fftw_plan_dft_2d(N, N, Phi, PhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_2d(N, N, NPhiF, NPhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_2d(N, N, Phi, Phi, FFTW_BACKWARD, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary); // read initial …

我用python编写了一个关于颜色检测的程序.但是'Erode'句子总是出现错误.这是我的计划的一部分.谢谢.

# Convert the image to a Numpy array since most cv2 functions
# require Numpy arrays.
frame = np.array(frame, dtype=np.uint8)

threshold = 0.05
#blur the image
frame=cv2.blur(frame, (5,5))
#Convert from BGR to HSV
hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
#split into 3
h, s, v= cv2.split(hsv)
#red color
s = cv2.threshold(h, 15, 1, cv2.THRESH_BINARY_INV)#1-15,x>15 y=0
h = cv2.threshold(h, 245, 1, cv2.THRESH_BINARY)#245-255 x>245 y=1
h = h + s
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3, 3)) 
h = cv2.erode(h, kernel)
v = cv2.sumElems(h)

我需要从二项式(n,p)分布生成随机数.

二项式(n,p)随机变量是n个均匀变量的总和,其中概率为1.在伪代码中,x=0; for(i=0; i<n; ++i) x+=(rand()<p?1:0);将生成二项式(n,p).

我需要为小的和非常大的n生成这个,例如n = 10 ^ 6和p = 0.02.是否有任何快速数值算法来生成它？

编辑 -

现在这是我的近似值(以及精确泊松和正态分布的函数) -

    public long Binomial(long n, double p) {
        // As of now it is an approximation
        if (n < 1000) {
            long result = 0;
            for (int i=0; i<n; ++i)
                if (random.NextDouble() < p) result++;
            return result;
        }
        if (n * p < 10) return Poisson(n * p);
        else if (n * (1 - p) < 10) return n - Poisson(n * …

var str = '0.25';

如何将上面的内容转换为 0.25?

我正在C#/ .NET中开展科学计算和可视化项目,我们使用doubles来表示所有物理量.由于浮点数总是包含一些舍入,我们有简单的方法来进行相等比较,例如:

static double EPSILON = 1e-6;

bool ApproxEquals(double d1, double d2) {
    return Math.Abs(d1 - d2) < EPSILON;
}

很标准.

然而,EPSILON当我们遇到"相等"数量的误差大于我们预期的情况时,我们不断发现自己必须调整幅度.例如,如果将5个大doubles 相乘然后除以5次,则会失去很多精度.它已经达到了我们无法使EPSILON过大的程度,否则它会给我们带来误报,但我们仍然会得到假阴性.

一般来说,我们的方法是寻找更多数值稳定的算法来使用,但程序是非常计算的,而且我们已经能够做到这一点.

有没有人有任何好的策略来处理这个问题？我Decimal稍微调查了这个类型,但我对性能感到担心,我对它的了解并不知道它是否可以解决问题或者只是模糊它.Decimal如果它能解决这些问题,我愿意接受一个适度的性能命中(比方说,2x),但性能肯定是一个问题,因为代码主要受到浮点运算的限制,我不认为它是一个无理的担忧.我看到人们引用了100倍的差异,这肯定是不可接受的.

此外,切换到Decimal其他复杂情况,例如Math库中普遍缺乏支持,因此我们必须编写自己的平方根函数,例如.

有什么建议？

编辑:顺便说一句,我使用常数epsilon(而不是相对比较)的事实不是我的问题.我只是把它作为一个例子,它实际上不是我的代码的snippit.改变相对比较不会对问题产生影响,因为问题是由于数字变得非常大而后来又变小而失去精确度.例如,我可能有一个值1000,然后我对它进行一系列计算,这将导致完全相同的数字,但由于精度损失,我实际上有1001.如果我然后去比较这些数字,它不会如果我使用相对或绝对比较(只要我以对问题和规模有意义的方式定义比较),那就更重要了.

无论如何,正如Mitch Wheat建议的那样,算法的重新排序确实有助于解决这些问题.

我必须从Python切换到C/C++.
你知道一个快速的"参考教程"或类似的东西,以引用如何开始？例如Numpy和Scipy教程.
例如,我已经阅读了很多"文档"

• C++ for dummies
• K&R C编程语言
• 许多博客和在线文档,例如:http://eli.thegreenplace.net/2010/01/11/pointers-to-arrays-in-c/,
• http://newdata.box.sk/bx/c/
• StackOverflow上的大量问答
• ...

但是我仍然不清楚如何开始移植到C/C++之类的东西:

#!/usr/bin/env python

import time
import numpy as np
import tables as tb

"""Retrieve 3D positions form 1000 files and store them in one single HDF5 file.
"""

t = time.time()

# Empty array
sample = np.array([])
sample.shape = (0,3)

# Loop over the files
for i in range(0, 1000):
  filename = "mill2sort-"+str(i)+"-extracted.h5"
  print "Doing ", filename
  # Open data file
  h5f = tb.openFile(filename, 'r')
  # Stack new …

我有一组点(x,y),我需要找到使用MATLAB通过原点的最佳拟合线.

当我在Python中使用小数字进行算术运算时,我得到以下意外结果:

>>> sys.float_info
sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)
>>> (1. - (1.e-17) ) < 1.
False

我知道浮点数不具有无限精度,但它应该能够处理像"1e-17"这样的"大"小数,不应该吗？