script_fmax.jl

#Julia 1.6.2
#RidePooling 0.8.2

using RidePooling
using DelimitedFiles

#error function (to be minimized)
function f(x,N,target,t0,specs)
    if x < 0
        return 1.0
    end

    model=get_model(;N_bus=N,ν=x/t0,specs...);

    #now comes the initial phase that will be skipped for the statistics
    run!(model,dropped=10*max(N,100),time=2*t0) #on average, every bus has dropped off >10 users + simulation time 4 times bigger than avrg. direct distance.
    startindex=length(model.requests)

    run!(model,requested=startindex+20*max(N,100))
    res=served_percentage(model,startindex=startindex)-target

    if abs(res)>0.5
        return res

    else
	run!(model,requested=startindex+50*max(N,100))
        res=served_percentage(model,startindex=startindex)-target
        if abs(res)>0.15
            return res

        else
	    run!(model,requested=startindex+100*max(N,100))
            res=served_percentage(model,startindex=startindex)-target
            if abs(res)>0.05
                return res

            else
		run!(model,requested=startindex+200*max(N,100))
                return served_percentage(model,startindex=startindex)-target
            end
        end
    end
end

function prnt(x)
    println(x)
    flush(stdout)
end

percentage=eval(Meta.parse(ARGS[1]))
index=eval(Meta.parse(ARGS[2]))

Ns=unique(Integer.(round.(exp.(range(log(1),log(1000);length=61)))))
N=Ns[index]
target=percentage/100

filepath=pwd()*"/../evaluation/f$(percentage)_$(N).txt"
guesspath=pwd()*"/../evaluation/guess_f$(percentage)_$(N).txt"

map_folder=pwd()*"/../../map/"
include(map_folder*"map.jl")
include(pwd()*"/../dispatcher.jl")

#initiate model
specs=(;
        map=map,
        route_matrix=RM,
        subspaces=:all_edges,
        routing=:lookup,
        speed_dict=speed_dict,
        seed=1
        )
specs=merge(specs,dispatcher)



TEST=false
init=1.6*N^1.14
occursin("32",pwd()) ? prnt("map with a 32. init prefactor = 1.6") : nothing
occursin("64",pwd()) ? (init*=2.6/1.6;prnt("map with a 64. init prefactor = 2.6")) : nothing
occursin("128",pwd()) ? (init*=5.5/1.6;prnt("map with a 128. init prefactor = 5.5")) : nothing


while TEST==false
    A=init*(0.95+0.1*rand())
    prnt("################")
    prnt("N = $(N)")
    prnt("initial guess A = $(round(A;digits=4))")
    prnt("################\n")
    fA=f(A,N,target,t0,specs)
    prnt("fA=$(round(fA;digits=4))")

    if fA>0
        prnt("A=$(round(A;digits=4)) too small, B must be larger.")
        best_positive=A
        best_negative=false
        B=A
        while best_negative==false
            B*=1.1
            prnt("Trying B=$(B)")
            fB=f(B,N,target,t0,specs)
            prnt("fB=$(fB)")
            best_negative= fB<0 ? B : false
            fB>0 ? ((A,fA)=(B,fB)) : nothing
        end

    else
        prnt("A too large, B must be smaller.")
        B,fB=(A,fA)
        prnt("Setting B=$(round(A;digits=4)) to maintain the convention that A<B")
        best_negative=B
        best_positive=false
        while best_positive==false
            A/=1.1
            prnt("Trying A=$(round(A;digits=4))")
            fA=f(A,N,target,t0,specs)
            prnt("fA=$(round(fA;digits=4))")
            best_positive= fA>0 ? A : false
            fA<0 ? ((B,fB)=(A,fA)) : nothing
            if A<1e-3
                prnt("This many buses don't seem to be able to serve $(percentage) percent of any demand. Exiting.")
                exit()
            end
        end
    end

    iteration=0
    while true
        if (fA<0.02 && fB>-0.02)
            prnt("")
            prnt("breaking free after $(iteration) iterations with")
            prnt("A=$(round(A;digits=4)), fA=$(round(fA;digits=4))")
            prnt("B=$(round(B;digits=4)), fB=$(round(fB;digits=4))")
            break
        end

        iteration+=1

        C_mean=(A+B)/2
        C_secante=B-fB*(B-A)/(fB-fA)
	    writedlm(guesspath,C_secante);prnt("\n  ..writing $(round(C_secante;digits=4)) (current secante root) to guess_path.")
        if fA<0.01
            C=(2*C_secante+B)/3
        elseif fB>0.01
            C=(2*C_secante+A)/3
        else
            C=(2*C_secante+C_mean)/3
        end
        fC=f(C,N,target,t0,specs)

                prnt("iteration $(iteration)")
                prnt("A=$(round(A;digits=4)), fA=$(round(fA;digits=4))")
                prnt("B=$(round(B;digits=4)), fB=$(round(fB;digits=4))")
                prnt("C=$(round(C;digits=4)), fC=$(round(fC;digits=4))")

        if fC>0
            A,fA=(C,fC)
        else
            B,fB=(C,fC)
        end
    end

    C=B-fB*(B-A)/(fB-fA)
	writedlm(guesspath,C);prnt("  ..writing $(round(C;digits=4)) (final secante root) to guess_path.")

    #find root of parabola
    fC=f(C,N,target,t0,specs)
    mat=[A^2 A 1;B^2 B 1;C^2 C 1]
    a,b,c=inv(mat)*[fA,fB,fC]
    p,q=(b/a,c/a)
    x1=-p/2-sqrt((p/2)^2-q)
    x2=-p/2+sqrt((p/2)^2-q)
    C= abs(x1-C)<abs(x2-C) ? x1 : x2
    writedlm(guesspath,C);prnt("  ..writing $(round(C;digits=4)) (root of parabola) to guess_path.")
                prnt("\nFINAL: $(round(C;digits=4))")

    #test served percentage
    fFINAL=f(C,N,target,t0,specs)

    global TEST=abs(fFINAL)<0.02
                prnt("TEST: served percentage = $(round(100*(fFINAL+target);digits=2))%\n\n\n")
    if TEST
      prnt("test worked.")
      writedlm(filepath,C)
    else
      global init=C
      prnt("test failed.")
    end
 end