//оптимизируем цену, стартуя из точки startPoint
//A,B,C - координаты городов
//prices - цены для городов A,B,C
long double Price(vector<long double> A, vector<long double> B, vector<long double> C, vector<long double> prices, vector<long double> startPoint)
{
vector<long double> center = //координаты фермы
startPoint, //точка начального приближения
grad = {1,1};
vector<long double> center1;
long double p1 = prices[0], p2 = prices[1], p3 = prices[2], //цены на строительство дорог из городов A,B,C
speed=0.0001, //скорость спуска
long double pCA = distance(center, A)*p1,
pCB = distance(center, B)*p2,
pCC = distance(center, C)*p3;
speed = 0.0001*std::max(pCA, std::max(pCB, pCC)); //делаем скорость соответствующей входным данным
price = 200000, newPrice = 200000, newPrice1 = 200000;
//цены, если поместить ферму в город
//если строить дорогу из этого города дорого или города на одной прямой, так бывает оптимально
long double prA = p2*sqrtl(powl(B[0] - A[0], 2) + powl(B[1] - A[1], 2)) +
p3*sqrtl(powl(C[0] - A[0], 2) + powl(C[1] - A[1], 2));
long double prB = p1*sqrtl(powl(A[0] - B[0], 2) + powl(A[1] - B[1], 2)) +
p3*sqrtl(powl(C[0] - B[0], 2) + powl(C[1] - B[1], 2));
long double prC = p1*sqrtl(powl(A[0] - C[0], 2) + powl(A[1] - C[1], 2)) +
p2*sqrtl(powl(B[0] - C[0], 2) + powl(B[1] - C[1], 2));
//минимальная цена, если поместить ферму в город
long double pr_point = std::min(prC, std::min(prA,prB));
long double sq1,sq2,sq3;
int iter = 0;
long double mul = 0;
int nm=0,nm1=1,nm2=2,iter_speed=0;
vector<vector<long double> > moves(3); //история 3 ходов
moves[0] = {0,0};
moves[1] = {1,11};
moves[2] = {223,2};
begin_:
do
{
price = newPrice;
sq1 = sqrtl(powl(A[0] - center[0], 2) + powl(A[1] - center[1], 2)); //расстояния от текущего положения фермы до городов
sq2 = sqrtl(powl(B[0] - center[0], 2) + powl(B[1] - center[1], 2));
sq3 = sqrtl(powl(C[0] - center[0], 2) + powl(C[1] - center[1], 2));
//пришли в город, тогда меняем скорость и начинаем сначала
if(sq1<=1e-15){
moves[0] = {1,8}; nm=0;nm1=1;nm2=2;moves[1] = {0,4}; moves[2] = {32,24}; iter = 0; newPrice = pr0; center = center0; speed = speed / 1.234234224; goto begin_;}
if(sq2<=1e-15){moves[0] = {1,8}; nm=0;nm1=1;nm2=2;moves[1] = {0,4}; moves[2] = {32,24}; iter = 0;newPrice = pr0;center = center0; speed = speed / 1.2624252523; goto begin_;}
if(sq3<=1e-15){moves[0] = {1,8}; nm=0;nm1=1;nm2=2;moves[1] = {0,4}; moves[2] = {32,24}; iter = 0;newPrice = pr0;center = center0; speed = speed / 1.2564333; goto begin_;}
long double grad01 = - p1*(A[0] - center[0]) / sq1;
long double grad02 = - p2*(B[0] - center[0]) / sq2;
long double grad03 = - p3*(C[0] - center[0]) / sq3;
long double grad11 = - p1*(A[1] - center[1]) / sq1;
long double grad12 = - p2*(B[1] - center[1]) / sq2;
long double grad13 = - p3*(C[1] - center[1]) / sq3;
grad[0] = grad01 + grad02 + grad03; //компоненты градиента цены по координатам фермы
grad[1] = grad11 + grad12 + grad13;
//пытаемся как можно больше продвинуться в направлении градиента
mul = 1000;
while(mul >= 1e-15)
{
center1 = center + grad * (-1.0) * mul * speed; //перемещаем ферму
long double sq11 = sqrtl(powl(A[0] - center1[0], 2) + powl(A[1] - center1[1], 2));
long double sq21 = sqrtl(powl(B[0] - center1[0], 2) + powl(B[1] - center1[1], 2));
long double sq31 = sqrtl(powl(C[0] - center1[0], 2) + powl(C[1] - center1[1], 2));
newPrice1 = p1 * sq11 + p2 * sq21 + p3 * sq31; //новая цена после перемещения фермы
if(newPrice1 < price)break;
if(mul == 0)break;
mul = mul / 2.;
}
center = center1;
newPrice = newPrice1;
pCA = distance(center, A)*p1; //цены которые получаются сейчас
pCB = distance(center, B)*p2;
pCC = distance(center, C)*p3;
++iter;
moves[(iter-1)%3] = center;
if(iter > 2)
{
nm = (iter - 1)%3;
if(nm<0)nm+=3;
nm1 = nm - 1;
if(nm1<0)nm1+=3;
nm2 = nm - 2;
if(nm2<0)nm2+=3;
//скалярное произведение ходов на шагах nm,nm-1,nm-2
if(dot(moves[nm1]-moves[nm2], moves[nm]-moves[nm1]) < 0) //перескочили через минимум туда и обратно
{
speed = speed / 2.; //уменьшаем скорость
}
}
++iter;
}
while(newPrice - price < 0); //двигаемся пока цена уменьшается
if(pr_point < price)return pr_point; //а вдруг лучше поместить ферму в один из городов?
return price;
}
//оптимизируем цену, стартуя из разных точек
long double optimizedPrice(vector<long double> A, vector<long double> B, vector<long double> C, vector<long double> prices)
{
long double P1 = Price(A,B,C,prices,(A+B)/2.); //середина сторона AB
long double P2 = Price(A,B,C,prices,(A+C)/2.); //середина сторона AC
long double P3 = Price(A,B,C,prices,(B+C)/2.); //середина сторона BC
long double P4 = Price(A,B,C,prices,(A+B+C)/3.); //центр треугольника
long double d = std::min(P4, std::min(P1,std::min(P2,P3)));
return d;
}