fix step size 过高引起的Precision异常

Projection method 里面分为两步，一步是步长测试，二是其他函数值得计算。从fix step size做起，可以先排除step size adjust 上的问题。
数据准备好，初次实验成功，下面看边界条件。
回想Fix 遇到什么问题？步长超过边界不能迭代
Ptt并不高
step size 没有下限，怎么低都可以，就是计算效率不高。比如极限情况，步长为0.0001，迭代数百万次，精度一直在提高。迭代数百万次，说明运行稳定。精度一直在提高说明，虽然慢，但没有错误。
下面看看另一个极端情况。

Step size = 10; Precision=e-5 Step size = 10; Precision=e-10 Step size = 10; Precision=e-15
发现Bug：精度突变。最后出现古怪符号。精度最高e-16
Step size = 10; Precision=e-20 Step size = 10; Precision=e-20

继续测试步长极限

Step size = 20; Precision=e-10 Step size = 30; Precision=e-10 Step size = 40; Precision=e-10

• Note: step size=40, Exception appears, the precision reduce and then, exception appears.

Step size = 50; Precision=e-10

• Note: step size=50, Exception again appears, the precision reduce and then, exception appears.

Step size = 60; Precision=e-10 Step size = 100; Precision=e-10 Step size = 1000; Precision=e-10

• note: we have not good results when step size or precision go to limitations.
  Route cost 没有问题
  Mode cost 在前面的测试中都是正常的， 这里也没有出现异常。只是数值比path cost要大一些。

Mode cost function

从程序上看不出问题

测试flow update，从结果上看不出问题，就是一个projection 操作，从中间过程看，也没有问题。首先是mode做流量调整，fix step size 乘以 transfer cost。投影。流量守恒。然后解路径流量。将升级后的流量分配到下辖的各条路径上。按mode demand plus 守恒。投影运算，流量-cost，没错。

然后计算误差

Paste_Image.png

误差没有缩小，反而变大了。但是结合以前的情况也没有问题。在接下来的几次迭代中，精度反复波动。第七次迭代，误差急剧变大。
第八次：not a number

path flow -- path cost (OK) -- mode cost (Exception) ---

检查mode cost
nest sum 出现有的nest一点流量都没分到的情况。
在上一步中所有流量都分到一个nest mode上。上一步就有问题了，流量分配不对。logit based 每种方式都有一定的流量。

Paste_Image.png 从上一步测试

ptt （OK）
mode cost ：
input是mode demand，其中一个demand是0，导致endogenous cost为负无穷。取值-999999999999999。外部cost是0.603289. 所以，总体cost是：

Paste_Image.png

然后，减去最小值就变成

Paste_Image.png

这个基本是没错的。有点误差。不知道为啥，小数点后面的问题。简单的减法，小数点后面几位不一致。

Paste_Image.png

正常情况下没有这个问题，可能是数值太大引起的。可以单独领出来试试。

总的来说就是Demand=0引起的。

把数值降下来就对了，减法就对了

那么就是怎么处理demand=0的情况？
在输入端控制，给一个极小的demand：e-20。

if (nm[nestmode_i].NestSum != 0) {
                nm[nestmode_i].EndogenousCost = \
                    (MuE / Theta)*log(0.00000000000000000001/ pow(nm[nestmode_i].Membership, 1 / MuE)) \
                    + ((1 - MuE) / Theta)*log(nm[nestmode_i].NestSum);
            }
            else {
                nm[nestmode_i].EndogenousCost = \
                    (MuE / Theta)*log(0.00000000000000000001 / pow(nm[nestmode_i].Membership, 1 / MuE)) \
                    + ((1 - MuE) / Theta)*log(0.00000000000000000002);
            }

Demand 调整以后，再测试精度。在第八次迭代出现not a number的情况。
在求ptt和path cost的时候应该没有问题了。
问题可能处在termination test上。
IndexNM[nestmode_i] = IndexNM[nestmode_i] + (p[path_i].Pf / nm[nestmode_i].Demand) * (p[path_i].PttTrans / p[path_i].Ptt);
分母nm[nestmode_i].Demand)确实可能为0；
至此，基本确定了fix step size 过高引起的Precision异常的问题。

Self-adaptive Projection 的问题，在ND上还是有问题，（Demand）

FW 解ND net 已经可以了。
Fix step size 至少处理了
S2FunValueNth(l, p, nm, od);
S3SolutionUpdate(p, nm, od, Alpha);
TerminationIndex = S8Termination(p, nm, od); printf("%.10lf\n", TerminationIndex);
S9Transit(p, nm);
都没什么问题，其中S2和S3是重要步骤，S8和S9相对简单。
在self adaptive中，还有更多的计算量。
2，3，4，5，6试算步长
7调整步长之Gamma调整
8收敛测试
9流量交替
9.1步长交替
先考虑那些东西要交替？flow/demand；alpha ；Gamma每次赋值就行。flow/demand cost都是计算而来的。在SAGP中涉及到的所有变量，x, y, alpha, alpha Plus, 这些变量中需要进行交替操作的只有input（流量）和步长。这两步骤很简单，且经过多次操作没有问题。
8收敛测试刚刚做过，没问题。
7步长调整中Gamma步长是步长调整阶段的变量。其他还有一个变量是最小整数l。S2-S6就是为了找到最优l，S7就是为了调整这个Gamma
S7和S6几乎相同，只是S6中的Delta是外部决定的，在S7直接给了0.5. 2-Delta始终大于0.5，所以S6中的不等式比S7更容易达到。如果S7中的不等式也成立，说明这个不等式条件很容易达成，步长可以扩大一点。由于S7和S6几乎相同，S7的检查也pass。
检查S6. 首先大于等于，基本公式都对过，没有问题。
S6中有三大函数：都是内积。内积的input分别是demand/flow；demand/flow cost；demand/flow cost （alpha）。先就这内积的形式来看看。三个内积的计算形式都很简单。没有问题，下面就是看看input对不对。
S5算demand/flow cost plus（alpha）
S4算demand/flow cost plus
S3算demand/flow plus
S2算demand/flow cost
S3 and S2 已经完成测试，没有太多问题。
S4 S5在一些情况下也正常工作，下面就极端情况做个别测试即可。FW就是在这种情况下测出问题的。

尝试跑起来

initial step size

when initial step size =5； step size does not change, but after 395 iterations, the precision collapse and go into endless loop

Paste_Image.png

when initial step size =10, the step size keeps increase, but after 94 iterations, the precision collapse and step size becomes zero, but didn't go into endless loop

Paste_Image.png

check above phenomenon

in 93th iteration, we obtain good results, and then these results will be like input.
the input is x, alpha k, gamma, the following y, alpha k plus,x plus, residual is the product of x, alpha k, and gamma.

input of nest iteration

so, the red part is the input variable, plus gamma and alpha k

1. we can see, Gamma=alpha k+1 / 0.9. match the definition.
   as definition, Part1 - Part2 >= Part3 holds.

   
   Part1 - Part2 - Part3
   

   to be more detail, it can be observed that part 1,2 and 3 are all small.

   
   part 1, part 2, part 3
2. input is mode demand/path flow, the output is path flow cost.
   path flow ---link flow ---link cost---path cost---- under each nest mode (path cost - min path cost)
3. input is min path cost under each mode, the output is nest mode cost
   min path cost under each mode + nest mode demand + dispersion parameter: Mu, Theta, Membership---nest mode cost
   Mode cost function
   the input is OK, nest mode demand is not equal to 0; so, the input is ok.
   demand, cost, endogenous demand, exogenous demand
   note: test whether we can use scientific notation for calculation. Yes, scientific notation can be used for calculation. --> mode cost, mode cost plus, we all use scientific notation to indicate a small enough value to replace zero.
4. solution update:

5. input: x, y, alpha k+1

   
   mode demad plus, mode demand, alpha*cost
   o,d,nest, mode,mode demad plus,mode demand, alpha*cost
   

   demand update is normal, flow conservation is held.

6. flow demand
7. y plus
   input is x plus,
   path flow plus---link flow ---link cost ---path cost plus---- under each nest mode (path cost plus - min path cost plus)
   min path cost under each mode plus + nest mode demand plus + dispersion parameter: Mu, Theta, Membership---nest mode cost plus
   目前为止都没问题。总结：Debug要回溯，逆序，从问题查起。

8. start from the problem.
   only one modification was made, give demand which supposed to be zero to be a very small number.
   in 95th iteration, the precision do not collapse, but the step size become zero.

   
   step size become a very small value.
   

   where is the problem? Input or operation?

9. see the operation in S6
   part 1 =0;
   part 2 ~~0;
   part 3 =0;
   before 95th iteration, part 1 !=0;
   so, check part 1.

   
   part 1, 2, 3
   part 1, 2, 3
   evaluation trend of part 1
   cost, cost P
   

   It is observed that cost p is unnormal. As expected, cost P in the previous iteration is stable.

   
   cost, cost P in previous iteration
   S4FunValueNPth(l, p, nm, od); produce unnormal cost p
   可能存在两个costP为0的情况
   确实是两个mode cost都为0，而mode cost P异常

为什么会有两个相同的最小值

Demand没有变异，正常，说明输入正常。
第一次接待就有问题，所以看第一次就好。

内层循环的第一次迭代
EndogenousCost
EndogenousCost + nm[nestmode_i].ExogenouCost 已经很接近了，如果算数运算正确的话，那确实相等。就是出现了两个最小值。跟我的理解，只能选一个。
如果第一个mode 的cost加上一个极小值，效果如何？
nm[0].Cost = nm[0].Cost + 1e-100; 1e-100太小，没有效果。
1e-15才有效果，1e-10太大进入死循环 也就是说精度一高，大小比价不出来。
这步就先不处理，其实本质上是精度处理能力不足引起的。

那解决的办法一般是用科学计数法来算了；

对于精度问题，想起前几天一件事，18个9减一个极小值，如0.15732，没有得出正确结果。 精度没有显示出来
15位9，确实有数字，但是结果不正确，小数点后面应该是8426

为什么cost P变异

第一次出现了流量不守恒的情况
输入没有问题，demandP，demand，alphak+1*TransCost
经过一次投影运算，flow确实没有问题

bug出现：存在多条最短路的时候，把非负最短路的流量加起来，就不对。只能选一条最短路。如果两条最短路都考虑，那么做流量升级的时候，没有好的升级方法来使得流量微调。所以只能选一条最短路。

解决办法，增加一个指标作为最短路指示指标。要做任何改动必须改一步试一步。

流量不守恒的问题得到纠正 然后又得到一个有趣的现象，精度是一个反复升降的过程。
度过多条最短路问题后，以后运行平稳，并且没有出现多条最短路问题，精度也持续提高

其实同理，在path flow update上也一样。

精度是一个反复升降的过程

精度的升降与是否存在多条最短路没有关系

从96与97次迭代看出，精度下降，但是并无多条最短路存在，说明精度下降的原因不在多条最短路上。
这里碰到的两条最短路也顺利过渡
这里的两条最短路也过渡顺利

总之，精度是一个反复升降的过程可能是在精度达到瓶颈以后遇到的问题，这里先不管。

精度最高到e-20

appendix

debug s2 s3

    for (path_i = 0; path_i < NoPs; path_i++) {
        //printf("%d %d %d %d %lf %lf\t%lf\t%lf\n", p[path_i].Path[0], p[path_i].Path[5], p[path_i].Nest, p[path_i].Mode, \
            p[path_i].Pf,  p[path_i].Ptt, p[path_i].PttTrans);
    }
    S22FunValueNth_ModeFun(nm, od);
    for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
        //printf("%d %d %d %d %lf %lf %lf %lf\n", nm[nestmode_i].O, nm[nestmode_i].D, \
            nm[nestmode_i].Nest, nm[nestmode_i].Mode, \
            nm[nestmode_i].Demand, nm[nestmode_i].Cost, nm[nestmode_i].EndogenousCost, nm[nestmode_i].ExogenouCost);
    }
    //S3SolutionUpdate(p, nm, od, Alpha);
    for (od_i = 0; od_i < NoODs; od_i++) {
        od[od_i].NMNonMinSum2 = 0;
        od[od_i].NMNonMinSum = 0;
        od[od_i].NMCounter = 0;
        od[od_i].Flag = 0;
    }
    for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
        if (nm[nestmode_i].CostTrans != 0) {
            nm[nestmode_i].DemandP = nm[nestmode_i].Demand - Alpha[1] * nm[nestmode_i].CostTrans;
            if (nm[nestmode_i].DemandP < 0) {
                nm[nestmode_i].DemandP = 0;
            }
            //printf("%d %d %d %d %.5e\t%.5e\t%.5e\n", nm[nestmode_i].O, nm[nestmode_i].D, nm[nestmode_i].Nest, nm[nestmode_i].Mode\
                , nm[nestmode_i].DemandP, nm[nestmode_i].Demand, Alpha[1] * nm[nestmode_i].CostTrans);
        }
    }
    for (od_i = 0; od_i < NoODs; od_i++) {
        for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
            if (nm[nestmode_i].CostTrans != 0 \
                && od[od_i].O == nm[nestmode_i].O && od[od_i].D == nm[nestmode_i].D) {
                od[od_i].NMNonMinSum = od[od_i].NMNonMinSum + nm[nestmode_i].DemandP;
            }
        }
    }
    for (od_i = 0; od_i < NoODs; od_i++) {
        for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
            if (nm[nestmode_i].CostTrans == 0 && od[od_i].NMCounter == 0 && \
                od[od_i].O == nm[nestmode_i].O && od[od_i].D == nm[nestmode_i].D) {
                nm[nestmode_i].DemandP = (od[od_i].Dem - od[od_i].NMNonMinSum);
                od[od_i].NMCounter = 1;
            }
        }
    }
    for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
        //printf("%d %d %d %d %.5e\n", nm[nestmode_i].O, nm[nestmode_i].D, nm[nestmode_i].Nest, nm[nestmode_i].Mode, nm[nestmode_i].DemandP);
    }