text:
module ScanForDuplicates
using CUDA, Logging,..CUDAGpuUtils,..WorkQueueUtils,..ScanForDuplicates, Logging,StaticArrays, ..IterationUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils
export @getIsToVal,@loadAndScanForDuplicates,@setIsToBeValidated, @scanForDuplicatesB,scanForDuplicatesMainPart,scanWhenDataInShmem,manageDuplicatedValue

# """
# after previous sync threads we already have the number of how much we increased number of results  relative to previous dilatation step
# now we need to go through  those numbers and in case some of the border queues were incremented we need to analyze those added entries to establish is there
# any duplicate in case there will be we need to decrement counter and set the corresponding duplicated entry to 0
# """
# macro scanForDuplicatesB(oldCount, countDiff,localOffset)#innerWarpNumb,shmemblockData,shmemSum,resListIndicies,metaData,xMeta,yMeta,zMeta,metaDataDims,localOffset,maxResListIndex,outerWarpLoop,alreadyCoveredInQueues,sourceShmem)
#     return esc(quote
#     #<13 becouse we are intrested only in border queues as only in those we can have overlap
#     if(innerWarpNumb<13)
#         sync_threads()
#         @exOnWarp innerWarpNumb begin
#            #we loaded data into threades registers so now we need to iteratively go through this and add this register info to shared memory to be available to all threads
#            @unroll for threadNumber in 1:32 # we need to analyze all thread id x
#                     if( threadIdxX()==threadNumber ) #now we need some  values that are in the registers  of the associated thread
#                         shmemSum[33,innerWarpNumb]=  $oldCount
#                         shmemSum[34,innerWarpNumb]=  $countDiff
#                         # #queue result offset
#                         shmemSum[35,innerWarpNumb]= $localOffset
#                         #will be used in order to keep track of proper new size of counter - will stay 0 if we have no duplicates
#                         shmemSum[36,innerWarpNumb]= 0
#                     end# if ( threadIdxX()==warpNumb )
#                 sync_warp()# now we should have the required number for scanning of new values for duplicates important that we are analyzing now given queue with just single warp - so we also get separate blocks as each block is represented by given idX
#                     scanForDuplicatesMainPart(shmemSum,innerWarpNumb,resListIndicies,metaData,xMeta,yMeta,zMeta,shmemblockData,metaDataDims,threadNumber,maxResListIndex,outerWarpLoop,alreadyCoveredInQueues, isInRange,resList,dataBdim)
#                 sync_threads()
#             end # for threadNumber
#         end #exOnWarp
#     end
# end)#qote
# end

# """
# main part of scanninf we are already on a correct warp; we already have old counter value and new counter value available in shared memory
# now we need to access the result queue starting from old counter
# """
# function  scanForDuplicatesMainPart(shmemSum,innerWarpNumb,resListIndicies,metaData,xMeta,yMeta,zMeta,shmemblockData,metaDataDims,threadNumber,maxResListIndex,outerWarpLoop,alreadyCoveredInQueues, isInRange,resList,dataBdim)
#     #this hold information wheather we have any new entries in a result queue
#     if(shmemSum[34,innerWarpNumb]>0 )
#         #as we can analyze 32 numbers at once if the amount of new results is bigger we need to do it in multiple passes
#         @unroll for scanIter in 0:fld(shmemSum[34,innerWarpNumb],32)
#                 # here we are loading data about linearized indicies of result in main  array depending on a queue we are analyzing it will tell about gold or other pas
#                     entryIndex = (shmemSum[33,innerWarpNumb]  +shmemSum[35,innerWarpNumb] + (scanIter*32) + threadIdxX())
#                     shmemSum[threadIdxX(),innerWarpNumb] = resListIndicies[entryIndex]
#                     #now becouse the entries from  paddings had not been saved into the dilatation arrays we need to do it now
#                     if(resList[7]==1) #checking wheather result entry originated from padding
#                         #becouse we are inside threadNumber loop we are sure that this datablock will not by modified by any other warp - the only thing we need to be sure is that two threads from the same warp will not try to modify the same integer ...
#                         #situation is diffrent for diffrent queues
#                             #in case of top and down we do not need to be afreaid of overwriting
#                             if(innerWarpNumb>8 && innerWarpNumb<13)

#                             end
#                             #for anterior posterior we need to be sure that xses are diffrent
#                             #for left and right we need to be sure that y are diffrent

#                         krowa

#                     end
#             sync_warp() # now we have 32 linear indicies loaded into the shared memory
#             #so we need to load some value into single value into thread and than go over all value in shared memory
#             scanWhenDataInShmem(shmemSum,innerWarpNumb, scanIter,resListIndicies,metaData,xMeta,yMeta,zMeta,metaDataDims,threadNumber,maxResListIndex,outerWarpLoop ,alreadyCoveredInQueues )
#             sync_warp()
#         end# for scanIter
#             ## all of those below happen only on one lane of each warp -  distributing work we nned to look for ideas to distribute work
#             #if this below it means that we have removed some entries from result list; we also check for thread to get appropriate metadata location
#             if(threadIdxX()==threadNumber)
#                 # the counbter reduced if needed of dupl;cated values
#                 updated = (shmemSum[33,innerWarpNumb]+shmemSum[34,innerWarpNumb])-shmemSum[36,innerWarpNumb]
#                 #alreadyCoveredInQueues will be later used to get global sums
#                 #added in this step
#                 if((shmemSum[34,innerWarpNumb]-shmemSum[36,innerWarpNumb]) >0)
#                     alreadyCoveredInQueues[innerWarpNumb]+=(shmemSum[34,innerWarpNumb]-shmemSum[36,innerWarpNumb])

#                 end
#                 #we set the counter to new value only if there were some subtractions done - some values were duplicated
#                 if((shmemSum[36,innerWarpNumb]>0 ) && isInRange )
#                     @setMeta(((getNewCountersBeg()) +innerWarpNumb) , updated)
#                   #krr metaData[xMeta+1,yMeta+1,zMeta+1,((getNewCountersBeg()) +innerWarpNumb) ]= updated

#                 end
#                 #it shoud NOT be analyzed if corrected counter value is the same as the total number for this queue
#                 # but we need to analyze this only if there was any new value added relatively to last pass
#                 if( (shmemSum[34,innerWarpNumb]-shmemSum[36,innerWarpNumb])  >0)
#                     # krr if(isInRange && updated!= metaData[(xMeta +1),(yMeta+1),(zMeta+1),getBeginingOfFpFNcounts()+innerWarpNumb ])

#                     if(isInRange && updated!= @accMeta((getBeginingOfFpFNcounts()+innerWarpNumb)) )
#                         #so in res shmem we have information weather we should  validate this queue ...
#                         shmemblockData[(threadIdxX())+(innerWarpNumb+21)*33]= 1
#                     end

#                  end
#             end
#              sync_warp()
#              #resetting the shmem values  that potentially could be not overwritten
#             for i in 1:3
#                 @ifX i if(shmemSum[36,innerWarpNumb]>0 )   shmemSum[33+i,innerWarpNumb]=0 end
#             end

#     end#if
# end

# """
# in this spot we already have 32 (not more at least ) values in shared memory
# we need now only to establish are there any duplicates
# """
# function scanWhenDataInShmem(shmemSum,innerWarpNumb, scanIter,resListIndicies,metaData,xMeta,yMeta,zMeta,metaDataDims,threadNumber,maxResListIndex ,outerWarpLoop ,alreadyCoveredInQueues)
#     @unroll for tempCount in (shmemSum[33,innerWarpNumb]):((shmemSum[33,innerWarpNumb]+shmemSum[34,innerWarpNumb])+1)
#         resListCurrIndex = tempCount +shmemSum[35,innerWarpNumb]
#        #now we need to make sure that we are not at spot whre this value is legitimite - so this is first occurence
#         if((resListCurrIndex!=( (shmemSum[33,innerWarpNumb]  +shmemSum[35,innerWarpNumb] + (scanIter*32) + threadIdxX()) ) ) && resListCurrIndex < maxResListIndex  )
#             #finally we iterate over all values in any given thread and compare to associated value in shared memory
#             if(shmemSum[threadIdxX(),innerWarpNumb]>0)
#                 #now we can make small optimazation and reduce global memory use and if reslist currindex is in range of currently loaded value in shared memory we will load them from shared memory
#                 shmemIndexBase = (shmemSum[33,innerWarpNumb]  +shmemSum[35,innerWarpNumb] + (scanIter*32))
#                 scannedVal = if(resListCurrIndex>shmemIndexBase && resListCurrIndex<(shmemIndexBase+32)  ) # range that is in the shmem
#                                 shmemSum[(resListCurrIndex- shmemIndexBase),innerWarpNumb]
#                             else
#                                 resListIndicies[resListCurrIndex]
#                             end

#                 #we need also to remember that we can have only 2 copies of the same result entry  we will keep only first one and second one we will remove
#                 if( ( scannedVal == shmemSum[threadIdxX(),innerWarpNumb]) &&  (tempCount>  ((scanIter*32) + threadIdxX()))   )
#                     #incrementing shared memory to later actualize counter
#                    CUDA.@atomic shmemSum[36,innerWarpNumb]+=1
#                     #if we have repeated value one entry in ids we set to 0
#                     resListIndicies[resListCurrIndex]=0
#                 end
#             end
#         end
#     end
# end #scanWhenDataInShmem

    """
    as we are operating under assumption that we do not know how many warps we have - we do not know the y dimension of thread block we need to load data into registers with a loop
    and within the same loop scan it for duplicates
    so if we have more than 12 warps we will execute the loop once - in case we have more we need to execute it in the loop
    iterThrougWarNumb - indicates how many times we need to  iterate to cover all 12 ques if we have at least 12 warps available (ussually we will) we will execute it once
    locArr,offsetIter,localOffset - variables used for storing some important constants

    """
    macro loadAndScanForDuplicates(iterThrougWarNumb,locArr,offsetIter)

        return esc(quote
            @unroll for outerWarpLoop in 0:$iterThrougWarNumb
                #represents the number of queue if we have enought warps at disposal it equals warp number so idY
                innerWarpNumb = (threadIdxY()+ outerWarpLoop*blockDimY())
                #at this point we have actual counters with correction for duplicated values  we can compare it with the  total values of fp or fn of a given queue  if we already covered
                #all points of intrest there is no point to futher analyze this block or padding

                @setIsToBeValidated()

                if(innerWarpNumb<15)
                    shmemSum[34,innerWarpNumb]= 0
                    shmemSum[35,innerWarpNumb]=0
                    shmemSum[36,innerWarpNumb]=0
                end
            end#for
            @exOnWarp 15 if(isInRange)

            @setMeta((getIsToBeAnalyzedNumb()+15), (@getIsToVal(1) || @getIsToVal(3)|| @getIsToVal(5)|| @getIsToVal(7)|| @getIsToVal(11)|| @getIsToVal(13)) )#sourceShmem[(threadIdxX())+33*8]
            end
            @exOnWarp 16 if(isInRange)

                @setMeta((getIsToBeAnalyzedNumb()+16 ),(@getIsToVal(2) || @getIsToVal(4)|| @getIsToVal(6)|| @getIsToVal(8)|| @getIsToVal(10)|| @getIsToVal(12)) )#sourceShmem[(threadIdxX())+33*6]
            end
    sync_threads()

end)#quote

        end #loadAndScanForDuplicates

        """
        loads value from res shmem about weather a queue with supplied numb has anything worth validating
        """
        macro getIsToVal(numb)
            return esc(quote
            (shmemPaddings[(threadIdxX())+($numb+21)*33])
        end)
        end

        """
        here we will mark in metadata weather there is anything to be verified - here in given que ie - weather it is possible in given queue to cover anything more in next dilatation step
        so it is important for analysisof this particular block is  there is true - is there is non 0 amount of points to cover in any queue of the block
        simultaneously the border queues should indicate for neighbouring blocks (given they exist ) is there is any point in analyzing the paddings ...
           so we need this data in 3 places
           1) for the getIsToBeAnalyzedNumb value in metadata of blocks around
           2) for isNotTobeAnalyzed for a current block
           3) also we should accumulate values pf counters - add and reduce across all blocks of tps and fps covered - so we will know when to finish the dilatation steps
     """
       macro setIsToBeValidated()
     return esc(quote
         @exOnWarp innerWarpNumb begin
             if(innerWarpNumb<15 && isInRange)
                #we need also to remember that data wheather there are any futher points of intrest is not only in the current block
                # so here we establish what are the coordinates of metadata of intrest so for example  our left fp and left FN are of intrest to block to the left ...
                newXmeta = xMeta+ (-1 * (innerWarpNumb==1 || innerWarpNumb==2)) + (innerWarpNumb==3 || innerWarpNumb==4)+1
                newYmeta = yMeta+ (-1 * (innerWarpNumb==5 || innerWarpNumb==6)) + (innerWarpNumb==7 || innerWarpNumb==8)+1
                newZmeta = zMeta+ (-1 * (innerWarpNumb==9 || innerWarpNumb==10)) + (innerWarpNumb==11 || innerWarpNumb==12)+1
                # #check are we in range
                if(newXmeta>0 && newYmeta>0 && newZmeta>0 && newXmeta<=metaDataDims[1] && newYmeta<=metaDataDims[2]  && newZmeta<=metaDataDims[3] && innerWarpNumb<13 )
                    metaData[newXmeta,newYmeta,newZmeta,getIsToBeAnalyzedNumb()+innerWarpNumb ] = (shmemPaddings[(threadIdxX())+(innerWarpNumb+21)*33] )

                end #if in meta data range
             end#if
             end#ex on warp

                #here we set the information weather any of the queue related to fp or fn in a particular block  has still something to be analyzed
     end)#quote
     end#setIsToBeValidated

end#module

# \
# """
# as we are operating under assumption that we do not know how many warps we have - we do not know the y dimension of thread block we need to load data into registers with a loop
# and within the same loop scan it for duplicates
# so if we have more than 12 warps we will execute the loop once - in case we have more we need to execute it in the loop
# iterThrougWarNumb - indicates how many times we need to  iterate to cover all 12 ques if we have at least 12 warps available (ussually we will) we will execute it once
# locArr,offsetIter,localOffset - variables used for storing some important constants

# """
# macro loadAndScanForDuplicates(iterThrougWarNumb,locArr,offsetIter)

#     return esc(quote

#     # @unroll for outerWarpLoop in 0:$iterThrougWarNumb
#     #         #represents the number of queue if we have enought warps at disposal it equals warp number so idY
#     #         innerWarpNumb = (threadIdxY()+ outerWarpLoop*blockDimY())
#     #         #now we will load the diffrence between old and current counter if we are in range of metadata
#     #         if(( innerWarpNumb)<15)
#     #             @exOnWarp innerWarpNumb begin
#     #                 if(isInRange)
#     #                     #store result in registers
#     #                     #old count
#     #                     $locArr = @accMeta((getOldCountersBeg()) +innerWarpNumb)
#     #                     # #diffrence new - old
#     #                     $offsetIter=   @accMeta(((getNewCountersBeg()) +innerWarpNumb) )- $locArr
#     #                     # #offset to find where are the results associated with given queue
#     #                     $localOffset =  @accMeta((getResOffsetsBeg()+innerWarpNumb))+$locArr
#     #                     # # enable access to information is it bigger than 0 to all threads in block
#     #                     shmemblockData[(threadIdxX())+(innerWarpNumb+6)*33] = UInt32($offsetIter >0)
#     #                     #for all border queues this information will be written down after scanning for duplicates here we are just getting info from main part ...
#     #                 else# if not in range
#     #                     $locArr = 0
#     #                     $offsetIter=0
#     #                     $localOffset=0
#     #                 end#if in range
#     #             end #@exOnWarp
#     #         end#if

#     #         sync_threads()
#     #         # so queaue 13 or 14
#     #             @exOnWarp innerWarpNumb begin
#     #                 if(innerWarpNumb==13 || innerWarpNumb==14)# so queaue 13 or 14
#     #                     if($offsetIter>0)
#     #                             # krr if($offsetIter!= metaData[xMeta+1,yMeta+1,zMeta+1,getBeginingOfFpFNcounts()+innerWarpNumb ])
#     #                         if($offsetIter!=  @accMeta((getBeginingOfFpFNcounts()+innerWarpNumb)) )
#     #                             shmemblockData[(threadIdxX())+(innerWarpNumb+21)*33]= 1
#     #                             #sourceShmem[(threadIdxX())+33*(6+(isodd(innerWarpNumb) *2) )]= true
#     #                         end
#     #                     end
#     #                     locOffset= UInt16(1)
#     #                     #setting what we need to locArr and reducing value in a warp
#     #                     $locArr+=$offsetIter
#     #                     #TODO (try to do warp reduction below instead of atomic... )
#     #                     if(($offsetIter)>0)
#     #                        CUDA.@atomic alreadyCoveredInQueues[innerWarpNumb]+=($offsetIter)
#     #                     end
#     #                     #@redOnlyStepOne(locOffset, shmemSum, $locArr, +)
#     #                     #now we have warp reduced value on first thread
#     #                     # @ifX 1 alreadyCoveredInQueues[innerWarpNumb]+=$locArr
#     #                 end#if
#     #             end#ex on warp
#     #             sync_threads()

#     #         #main function for scanning

#     #        @scanForDuplicatesB($locArr, $offsetIter,$localOffset)#$locArr, $offsetIter,innerWarpNumb,shmemblockData,shmemSum,resListIndicies,metaData,xMeta,yMeta,zMeta,metaDataDims,$localOffset,maxResListIndex,outerWarpLoop,alreadyCoveredInQueues,sourceShmem)

#     #        if(innerWarpNumb<15)
#     #             shmemSum[threadIdxX(),innerWarpNumb]=0
#     #         end

#     #     end#for
#     #     sync_threads()

#         @unroll for outerWarpLoop in 0:$iterThrougWarNumb
#             #represents the number of queue if we have enought warps at disposal it equals warp number so idY
#             innerWarpNumb = (threadIdxY()+ outerWarpLoop*blockDimY())
#             #at this point we have actual counters with correction for duplicated values  we can compare it with the  total values of fp or fn of a given queue  if we already covered
#             #all points of intrest there is no point to futher analyze this block or padding

#             @setIsToBeValidated()

#             if(innerWarpNumb<15)
#                 shmemSum[34,innerWarpNumb]= 0
#                 shmemSum[35,innerWarpNumb]=0
#                 shmemSum[36,innerWarpNumb]=0
#             end
#         end#for
#         @exOnWarp 15 if(isInRange)

#         @setMeta((getIsToBeAnalyzedNumb()+15), (@getIsToVal(1) || @getIsToVal(3)|| @getIsToVal(5)|| @getIsToVal(7)|| @getIsToVal(11)|| @getIsToVal(13)) )#sourceShmem[(threadIdxX())+33*8]
#         end
#         @exOnWarp 16 if(isInRange)

#             @setMeta((getIsToBeAnalyzedNumb()+16 ),(@getIsToVal(2) || @getIsToVal(4)|| @getIsToVal(6)|| @getIsToVal(8)|| @getIsToVal(10)|| @getIsToVal(12)) )#sourceShmem[(threadIdxX())+33*6]
#         end
# sync_threads()

# end)#quote


text:
@testset "Magnetostatics Functional Tests" begin
    @testset "Repelling magnetic dipoles" begin
        d1 = 0.01 # m
        m1 = 5e-6 # kg
        m2 = 5e-6 # kg
        ρ = 7800 # kg/m^3
        M = 1.2e6 # A/m
        mm1 = SVector(0.0, 0.0, M * m1 / ρ)
        mm2 = SVector(0.0, 0.0, M * m2 / ρ)
        μ_4π = 1e-7

        t1 = 0.0  # s
        t2 = 1.0 # s
        τ = (t2 - t1) / 100

        p1 = MagneticParticle(SVector(-d1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m1, mm1)
        p2 = MagneticParticle(SVector(d1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m2, mm2)
        parameters = MagnetostaticParameters(μ_4π)
        system = PotentialNBodySystem([p1, p2], Dict(:magnetic => parameters))
        simulation = NBodySimulation(system, (t1, t2))
        sim_result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        d2 = norm(get_position(sim_result, t2, 2) - get_position(sim_result, t2, 1))
        v_expected = sqrt(μ_4π * ( dot(mm1, mm2) * (1 / d1^3 - 1 / d2^3) ) / m1)
        v_actual = norm(get_velocity(sim_result, t2, 2))

        ε = 0.001
        @test v_expected ≈ v_actual atol = ε
    end

    @testset "Constructing magnetostatic potential parameters entity" begin
        default_potential = MagnetostaticParameters()
        @test 1e-7 == default_potential.μ_4π

        io = IOBuffer()
        potential1 = MagnetostaticParameters()
        potential2 = MagnetostaticParameters(4π)
        @test sprint(io -> show(io, potential1)) == "Magnetostatic:\n\tμ/4π:1.0e-7\n"
        @test sprint(io -> show(io, potential2)) == "Magnetostatic:\n\tμ/4π:12.566370614359172\n"
    end
end

ext: jl


text:
function integral(A::AbstractInterpolation, t::Number)
    bw, fw = samples(A)
    idx = max(1+bw, min(searchsortedlast(A.t, t), length(A.t) - fw))
    _integral(A, idx, t)
end

function integral(A::AbstractInterpolation, t1::Number, t2::Number)
    bw, fw = samples(A)
    # the index less than or equal to t1
    idx1 = max(1+bw, min(searchsortedlast(A.t, t1), length(A.t) - fw))
    # the index less than t2
    idx2 = max(2+bw, min(searchsortedlast(A.t, t2), length(A.t) - fw))
    if A.t[idx2] == t2
        idx2-=1
    end
    total = zero(eltype(A))
    for idx in idx1:idx2
        lt1 = idx == idx1 ? t1 : A.t[idx]
        lt2 = idx == idx2 ? t2 : A.t[idx+1]
        total += _integral(A, idx, lt2)-_integral(A, idx, lt1)
    end
    total
end

samples(A::LinearInterpolation{<:AbstractVector}) = (0, 1)
function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number)
  t1 = A.t[idx]
  t2 = A.t[idx+1]
  u1 = A.u[idx]
  u2 = A.u[idx+1]
  t^2*(u1 - u2)/(2*t1 - 2*t2) + t*(t1*u2 - t2*u1)/(t1 - t2)
end

samples(A::ConstantInterpolation{<:AbstractVector}) = (0, 1)
function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
  if A.dir === :left
    # :left means that value to the left is used for interpolation
    return A.u[idx]*t
  else
    # :right means that value to the right is used for interpolation
    return A.u[idx+1]*t
  end
end

samples(A::QuadraticInterpolation{<:AbstractVector}) = (0, 2)
function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number)
  t1 = A.t[idx]
  t2 = A.t[idx+1]
  t3 = A.t[idx+2]
  u1 = A.u[idx]
  u2 = A.u[idx+1]
  u3 = A.u[idx+2]
  (t^3*(-t1*u2 + t1*u3 + t2*u1 - t2*u3 - t3*u1 + t3*u2)/
  (3*t1^2*t2 - 3*t1^2*t3 - 3*t1*t2^2 + 3*t1*t3^2 + 3*t2^2*t3 - 3*t2*t3^2) +
  t^2*(t1^2*u2 - t1^2*u3 - t2^2*u1 + t2^2*u3 + t3^2*u1 - t3^2*u2)/
  (2*t1^2*t2 - 2*t1^2*t3 - 2*t1*t2^2 + 2*t1*t3^2 + 2*t2^2*t3 - 2*t2*t3^2) +
  t*(t1^2*t2*u3 - t1^2*t3*u2 - t1*t2^2*u3 + t1*t3^2*u2 + t2^2*t3*u1 - t2*t3^2*u1)/
  (t1^2*t2 - t1^2*t3 - t1*t2^2 + t1*t3^2 + t2^2*t3 - t2*t3^2))
end

samples(A::QuadraticSpline{<:AbstractVector{<:Number}}) = (0, 1)
function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
  t1 = A.t[idx]
  t2 = A.t[idx+1]
  u1 = A.u[idx]
  z1 = A.z[idx]
  z2 = A.z[idx+1]
  t^3*(z1 - z2)/(6*t1 - 6*t2) + t^2*(t1*z2 - t2*z1)/(2*t1 - 2*t2) + t*(-t1^2*z1 - t1^2*z2 + 2*t1*t2*z1 + 2*t1*u1 - 2*t2*u1)/(2*t1 - 2*t2)
end

samples(A::CubicSpline{<:AbstractVector{<:Number}}) = (0, 1)
function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
  t1 = A.t[idx]
  t2 = A.t[idx+1]
  u1 = A.u[idx]
  u2 = A.u[idx+1]
  z1 = A.z[idx]
  z2 = A.z[idx+1]
  h2 = A.h[idx+1]
  (t^4*(-z1 + z2)/(24*h2) + t^3*(-t1*z2 + t2*z1)/(6*h2) +
  t^2*(h2^2*z1 - h2^2*z2 + 3*t1^2*z2 - 3*t2^2*z1 - 6*u1 + 6*u2)/(12*h2) +
  t*(h2^2*t1*z2 - h2^2*t2*z1 - t1^3*z2 - 6*t1*u2 + t2^3*z1 + 6*t2*u1)/(6*h2))
end


text:
## CD3 γϵ and δϵ cooperative binding to TCRαβ ecotodomain
##

# bimolecular dissociations
#       γϵ + TCR ← γϵ-TCR
#       δϵ + TCR ← δϵ-TCR

##
using Revise
include("model_fit_global.jl")
include("models.jl")
include("helper_functions.jl")
include("binning_helper_fs.jl")
using DataFrames, Plots
CD3_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/dissociations.csv",DataFrame)

γϵ = delete_missing(CD3_data[:,2:3])
δϵ = delete_missing(CD3_data[:,8:9])
mix = delete_missing(CD3_data[:,14:15])
t_γϵ_thrml = delete_missing(CD3_data[:,5:5])
t_δϵ_thrml = delete_missing(CD3_data[:,11:11])
t_mix_thrml = delete_missing(CD3_data[:,17:17])
p_γϵ_thrml = survival_p.(t_γϵ_thrml)
p_δϵ_thrml = survival_p.(t_δϵ_thrml)
p_mix_thrml = survival_p.(t_mix_thrml)

s1=scatter(γϵ[1],γϵ[2], label="γϵ")
s2=scatter(δϵ[1],δϵ[2], label="δϵ")
s3=scatter(mix[1],mix[2], label="mix")
plot(s1,s2,s3)
plotly()
plot(t_γϵ_thrml,normalize_lnp.(p2lnp(p_γϵ_thrml)), label="γϵ")
plot!(t_δϵ_thrml,normalize_lnp.(p2lnp(p_δϵ_thrml)), label="δϵ")
plot!(t_mix_thrml,normalize_lnp.(p2lnp(p_mix_thrml)), label="mix")

## Thermal fluctuation bimolecular dissociation

tspan = (0.0,10.0)
t = collect(0.0:0.1:10)

gr()

u₀ = [1.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
# model_opt = (-5.,3.)
model_opt = (0.,100.)
pinit =[0.0]
dens = fill([],length(t_γϵ_thrml))
f_thrml = [[0.0]]
model = bimolecular_diss!

# CD3γϵ ------------------------------------------------------------------------
rates_γϵ_thrml,hes_γϵ_thrml,l_γϵ_thrml,boot_γϵ_thrml = optimization(model,p_γϵ_thrml,
                    dens,f_thrml,model_opt,t_γϵ_thrml)

# p = 10 .^ rates_γϵ_thrml
p = rates_γϵ_thrml
# p = [14.1]
k₋₁ᵒ = rates_γϵ_thrml
show_diss_fit(model,t_γϵ_thrml,p2lnp(p_γϵ_thrml),f_thrml)

# CD3δϵ ------------------------------------------------------------------------
rates_δϵ_thrml,hes_δϵ_thrml,l_δϵ_thrml,boot_δϵ_thrml = optimization(model,p_δϵ_thrml,
                    dens,f_thrml,model_opt,t_δϵ_thrml)
p = rates_δϵ_thrml
# p = [11.1]
k₋₂ᵒ = rates_δϵ_thrml
show_diss_fit(model,t_δϵ_thrml,p2lnp(p_δϵ_thrml),f_thrml)

# mix ------------------------------------------------------------------------
cons_mix_thrml = [[k₋₁ᵒ[1],k₋₂ᵒ[1]]]
cons_mix_thrml = [[14.1,11.1]]
u₀ = [0.4,0.4,0.2]          # Initial condition
model = cd3_dissociation!
rates_mix_thrml,hes_mix_thrml,l_mix_thrml,boot_mix_thrml = optimization(model,p_mix_thrml,
                    dens,cons_mix_thrml,model_opt,t_mix_thrml)
p = rates_mix_thrml
k₋₃ᵒ = rates_δϵ_thrml[1]
# f = (du,u,p,t) -> cd3_dissociation!(du,u,p,t,[],cons_mix_thrml[1])
# sol = solve_model(k₋₃ᵒ)
# sumsol =  sum(sol,dims=1)
show_diss_fit(model,t_mix_thrml,p2lnp(p_mix_thrml),cons_mix_thrml)

## Load Force-dependent bimolecular dissociation

f_γϵ, t_γϵ, γϵ_bin_ids = bin_dissociations(γϵ,6;bin_type="force")
f_δϵ, t_δϵ, δϵ_bin_ids = bin_dissociations(δϵ,6)
f_mix, t_mix, mix_bin_ids = bin_dissociations(mix,6)
f_γϵ, t_γϵ = resort(f_γϵ, t_γϵ)
f_δϵ, t_δϵ = resort(f_δϵ, t_δϵ)
f_mix, t_mix = resort(f_mix, t_mix)
p_γϵ = survival_p.(t_γϵ)
p_δϵ = survival_p.(t_δϵ)
p_mix = survival_p.(t_mix)

s1=scatter(f_γϵ, t_γϵ, title="γϵ")
s2=scatter(f_δϵ, t_δϵ, title="δϵ")
s3=scatter(f_mix, t_mix, title="mix")
plot(s1,s2,s3)

s1=scatter(t_γϵ, p2lnp(p_γϵ), title="γϵ")
s2=scatter(t_δϵ, p2lnp(p_δϵ), title="δϵ")
s3=scatter(t_mix, p2lnp(p_mix), title="mix")
plot(s1,s2,s3)

using StatsBase
t_γϵ_avg = mean.(t_γϵ)
f_γϵ_avg = mean.(f_γϵ)
plot(f_γϵ_avg,t_γϵ_avg)

t_γϵ[1]
prob = survival_p_2(t_γϵ[1])

survival_p_2(t) = collect((length(t)-1 :-1:0) / length(t))
plot(t_γϵ[1],p2lnp(prob))
mean(t_γϵ[2])
sum(t_γϵ[1] .* prob ./ sum(prob))
sum(t_γϵ[2].*survival_p_2(t_γϵ[2]))

# u₀ = [1.0]          # Initial condition
# tspan = (0.0,15.0)  # Simulation interval
# # model_opt = (-5.,3.)
# model_opt = (0.,100.)
# pinit =[1.0,1.0]
# dens = fill([],length(t_γϵ))
# model = bimolecular_diss_bell_single!
# cons_γϵ = [[k₋₁ᵒ[1], f] for f in f_γϵ_avg]
#
# rates_γϵ,hes_γϵ,l_γϵ,boot_γϵ = optimization(model,p_γϵ,
#                     dens,f_γϵ_avg,model_opt,t_γϵ)
#
# p = rates_γϵ
# show_diss_fit(model, t_γϵ, p2lnp(p_γϵ), f_γϵ_avg)
# plot_F_vs_t(f_γϵ,t_γϵ)

## try bimolecular dissociation bell model with activated state

u₀ = [1.0,0.0]
tspan = (0.0,15.0)
# model_opt = (-5.,3.)  # Simulation interval
model_opt = (0.,100.)
pinit = ones(8)
dens = fill([],length(t_γϵ))
model = bimolecular_diss_bell!
rates_γϵ,hes_γϵ,l_γϵ,boot_γϵ = optimization(model,p_γϵ,
                    dens,f_γϵ_avg,model_opt,t_γϵ)

p = rates_γϵ
# p = 10 .^ rates_γϵ
show_diss_fit(model, t_γϵ, p2lnp(p_γϵ), f_γϵ_avg)
Fs = collect(0.:0.5:30.)
plot()
plot_F_vs_t(f_γϵ,t_γϵ)
plot_force_lifetime_fit(model,Fs,f_γϵ,t_γϵ)

## try to include thermal fluctuation data to bfp data

p_γϵ_all = vcat(p_γϵ_thrml, p_γϵ)
t_γϵ_all = vcat(t_γϵ_thrml, t_γϵ)
f_γϵ_thrml = [zeros(size(t_γϵ_thrml[1]))]
f_γϵ_all = vcat(f_γϵ_thrml,f_γϵ)
f_γϵ_avg_all = vcat(0.0,f_γϵ_avg)
dens = fill([],length(t_γϵ_all))

rates_γϵ,hes_γϵ,l_γϵ,boot_γϵ = optimization(model,p_γϵ_all,
                    dens,f_γϵ_avg_all,model_opt,t_γϵ_all)

p = rates_γϵ
# p = 10 .^ rates_γϵ
show_diss_fit(model, t_γϵ_all, p2lnp(p_γϵ_all), f_γϵ_avg_all)
Fs = collect(0.:0.5:30.)
plot()
plot_F_vs_t(f_γϵ_all,t_γϵ_all)
plot_force_lifetime_fit(model,Fs,f_γϵ,t_γϵ)


text:
#=

Filename: Model.jl.
Author: Shunsuke Hori

Defining an instance "Model" storing parameters, functional forms, and grid.

=#

type Model
    alpha::Float64    # Productivity parameter
    beta::Float64     # Discount factor
    gamma::Float64    # risk aversion
    mu::Float64       # First parameter in lognorm(mu, sigma)
    s::Float64        # Second parameter in lognorm(mu, sigma)
    grid_min::Float64 # Smallest grid point
    grid_max::Float64 # Largest grid point
    grid_size::Signed # Number of grid points
    u::Function       # utility function
    u_prime::Function # derivative of utility function
    f::Function       # production function
    f_prime::Function # derivative of production function
    grid::Vector{Float64}   # grid
end

"""
construct Model instance using the information of parameters and functional form

arguments: see above

return: Model type instance
"""
function Model(; alpha::Float64=0.65,   # Productivity parameter
                 beta::Float64=0.95,    # Discount factor
                 gamma::Float64=1.0,    # risk aversion
                 mu::Float64=0.0,       # First parameter in lognorm(mu, sigma)
                 s::Float64=0.1,        # Second parameter in lognorm(mu, sigma)
                 grid_min::Float64=1e-6,# Smallest grid point
                 grid_max::Float64=4.0, # Largest grid point
                 grid_size::Signed=200, # Number of grid points
                 u::Function= c->(c^(1-gamma)-1)/(1-gamma), # utility function
                 u_prime::Function = c-> c^(-gamma),        # u'
                 f::Function = k-> k^alpha,                 # production function
                 f_prime::Function = k -> alpha*k^(alpha-1) # f'
                 )

    grid=collect(linspace(grid_min, grid_max, grid_size))

    if gamma == 1   # when gamma==1, log utility is assigned
        u_log(c) = log(c)
        m = Model(alpha, beta, gamma, mu, s, grid_min, grid_max,
                grid_size, u_log, u_prime, f, f_prime, grid)
    else
        m = Model(alpha, beta, gamma, mu, s, grid_min, grid_max,
                  grid_size, u, u_prime, f, f_prime, grid)
    end
    return m
end


text:
using GeoPhyInv

mgrid=repeat(range(-1000.0,stop=1000.0,length=201), 2) # test spatial grid
#ageom_surf=GeoPhyInv.Gallery.AGeom(mgrid,:surf,nss=5,nr=30);


text:
include(joinpath(@__DIR__, "..", "src", "artifact_funcs.jl"))

function trigger_download()
    output_artifact_path = get_path_to_artifact()
    return nothing
end
trigger_download()


text:
module MuttsTest

using Mutts
using Test
using MacroTools

#mutable struct Foo <: Mutt
#    x :: Int
#    y :: Int
#    __mutt_mutable :: Bool
#end
#
#Foo(x :: Int, y :: Int) = Foo(x, y, true)
#Foo(f :: Foo) = Foo(f.x, f.y)
#
#f = Foo(3,5)
#@assert ismutable(f)
#f.x = 4
#
#Base.copy(v :: Foo) = Foo(v.x, v.y)
#g = branch!(f)
#@assert !ismutable(f)
#@assert ismutable(g)
#markimmutable!(g)
#
#@assert !ismutable(g)

#mutable struct Bar <: Mutt
#    f :: Foo
#    z :: Float64
#    __mutt_mutable :: Bool
#end
#
#Bar(f :: Foo, z :: Float64) = Bar(f, z, true)
#Bar(b :: Bar) = Bar(f, z)
#
#

# Mutts macro expansion

@testset "macro parsing" begin
    @test @macroexpand(@mutt struct S end) isa Expr
    @test @macroexpand(@mutt struct S x end) isa Expr

    # Bad parses
    @test_throws Exception @macroexpand(@mutt 1)
    @test_throws Exception @macroexpand(@mutt begin
        struct A
        end
        struct B
        end
    end)
end
@testset "Simple Macro Usage" begin
    @eval begin
        @mutt struct S
            x
        end
        s = S(1)
        @test ismutable(s)
        markimmutable!(s)
        @test !ismutable(s)

        Base.copy(s::S) = S(s.x)

        s = S(1)
        @test ismutable(branch!(s))
        @test !ismutable(s)

        # Assignment from markimmutable! (PR #4)
        s1 = markimmutable!(s)
        @test s1 == s
    end
end

@mutt struct M
    x
end
@testset "Can't modify immutable Mutts" begin
    m = M(1)
    m.x = 2
    @test m.x == 2
    markimmutable!(m)
    @test_throws Exception m.x = 3
end

@testset "Inner constructors" begin
    # (eval required to create structs inside a Testset)
    @eval begin
        # Test Inner Constructors
        @mutt struct SimpleFields
            x
            y
        end
        @mutt struct NoCustomInner
            x :: Int
            y
        end
        @mutt struct WithCustomInners
            x :: Int
            y
            WithCustomInners(x, y=2) = new(x,y)
        end
        function make end
        @mutt struct WithInnerFunctions
            x :: Int
            y
            WithInnerFunctions(x, y=2) = new(x,y)
            # keyword arguments constructor
            WithInnerFunctions(; x=1, y=2) = new(x,y)
            # Custom inner function
            @__MODULE__().make() = new(1,2)
        end
        @mutt struct DefaultsWithTypeParams{A,B}
            x :: A
            y :: B
        end
        @mutt struct CustomWithTypeParams{A,B}
            x :: A
            y :: B
            CustomWithTypeParams() = new{Int,Int}(1,2)
            CustomWithTypeParams{A,B}(x,y) where {A,B} = new{A,B}(x,y)
        end
    end
    @eval begin
        @test ismutable(SimpleFields(1,2))
        @test ismutable(NoCustomInner(1,2))
        @test ismutable(WithCustomInners(1,2))
        @test ismutable(WithCustomInners(1))

        @test ismutable(WithInnerFunctions(1))
        @test ismutable(WithInnerFunctions(y=1))
        @test ismutable(make())

        @test ismutable(DefaultsWithTypeParams{Int,String}(1,""))
        @test ismutable(CustomWithTypeParams())
        @test ismutable(CustomWithTypeParams{Int,String}(1,""))
    end
end

@testset "structs with supertypes" begin
    @eval begin
        @mutt struct MuttsFloat <: AbstractFloat
            v :: Float64
        end
        @test MuttsFloat <: AbstractFloat
        v = MuttsFloat(2.0)
        @test v isa AbstractFloat

        v.v = 3.0
        @test v.v == 3.0
        markimmutable!(v)
        @test_throws Exception v.v = 4.0
    end
end

end


text:
ExodusError = Int64

function exodus_error_check(error::ExodusError, method_name::String)
    if error < 0
        error("Error from exodus library call in method $method_name")
    end
end


text:
# Copyright 2015, 2016, 2017, 2018 Martin Holters
# See accompanying license file.

include("checklic.jl")

using ACME
using Compat: argmin, range
import Compat.Test
using Compat.Test: @test, @test_broken, @test_throws, @test_warn, @testset
using ProgressMeter
using Compat.SparseArrays: sparse, spzeros

@testset "topomat" begin
    tv, ti = ACME.topomat(sparse([1 -1 1; -1 1 -1]))
    @test tv*ti'==spzeros(2,1)

    # Pathological cases for topomat:
    # two nodes, one loop branch (short-circuited) -> voltage==0, current arbitrary
    @test ACME.topomat(spzeros(Int, 2, 1)) == (hcat([1]), spzeros(0, 1))
    # two nodes, one branch between them -> voltage arbitrary, current==0
    @test ACME.topomat(sparse([1,2], [1,1], [1,-1])) == (spzeros(0, 1), hcat([1]))
end

@testset "LinearSolver" begin
    solver = ACME.LinearSolver(3)
    A = [1.0 0.5 0.4; 2.0 4.0 1.7; 4.0 7.0 9.1]
    @test ACME.setlhs!(solver, A)
    x = rand(3)
    y = similar(x)
    ACME.solve!(solver, y, x)
    @test A*y ≈ x
    y = copy(x)
    ACME.solve!(solver, y, y)
    @test A*y ≈ x
    @test_throws DimensionMismatch ACME.setlhs!(solver, zeros(2, 3))
    @test_throws DimensionMismatch ACME.setlhs!(solver, zeros(3, 4))
    @test_throws DimensionMismatch ACME.setlhs!(solver, zeros(4, 4))
    @test_throws DimensionMismatch ACME.solve!(solver, zeros(2), zeros(3))
    @test_throws DimensionMismatch ACME.solve!(solver, zeros(3), zeros(4))
    @test_throws DimensionMismatch ACME.solve!(solver, zeros(4), zeros(4))
    @test !ACME.setlhs!(solver, zeros(3,3))
end

@testset "simple circuits" begin
    @testset "empty circuit" begin
        circ = @circuit begin end
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 20)) == zeros(0, 20)
    end

    @testset "only one resistor" begin
        circ = @circuit begin
            r = resistor(0), [1] ⟷ [2]
        end
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 20)) == zeros(0, 20)
    end

    # simple circuit: resistor and diode in series, driven by constant voltage,
    # chosen such that a prescribe current flows
    @testset "resistor-diode" begin
        i = 1e-3
        r=10e3
        is=1e-12
        v_r = i*r
        v_d = 25e-3 * log(i/is+1)
        circ = @circuit begin
            vsrc = voltagesource(v_r + v_d), [+] ⟷ "supply voltage",[-] ⟷ gnd
            r1 = resistor(r)
            d = diode(is=is), [-] ⟷ gnd, [+] ⟷ r1[2]
            vprobe = voltageprobe(), [-] ⟷ gnd, [+] ⟷ r1[2]
            r1[1] ⟷ "supply voltage"
        end
        model = DiscreteModel(circ, 1)
        y = run!(model, zeros(0, 1))
        @test y[1] ≈ v_d
    end
end

@testset "circuit manipulation" begin
    @testset "programmatic reconnection" begin
        circ = @circuit begin
            r1 = resistor(10)
            r2 = resistor(100), [1] ⟷ r1[1], [2] ⟷ r1[2]
            src = voltagesource(1), [-] ⟷ r1[2]
            probe = currentprobe(), [+] ⟷ src[+],
                                    [-] ⟷ r1[1]
        end
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/10 + 1/100
        disconnect!(circ, (:r2, 1))
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/10
        disconnect!(circ, (:r1, 2))
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 0
        connect!(circ, (:r1, 2), (:r2, 1))
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/(10+100)
    end

    @testset "element deletion" begin
        circ = @circuit begin
            r1 = resistor(10)
        end
        r2_des = add!(circ, resistor(100))
        add!(circ, :r3, resistor(470))
        r4_des = add!(circ, resistor(1000))
        add!(circ, :src, voltagesource(1))
        add!(circ, :probe, currentprobe())
        connect!(circ, (:src, :+), (:probe, :+))
        connect!(circ, (:probe, :-), (:r1, 1), (r2_des, 1), (:r3, 1), (r4_des, 1))
        connect!(circ, (:src, :-), (:r1, 2), (r2_des, 2), (:r3, 2), (r4_des, 2))
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/10 + 1/100 + 1/470 + 1/1000
        delete!(circ, :r1)
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/100 + 1/470 + 1/1000
        delete!(circ, r4_des)
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/100 + 1/470
        delete!(circ, :r3)
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 1/100
        delete!(circ, r2_des)
        model = DiscreteModel(circ, 1)
        @test run!(model, zeros(0, 1))[1,1] ≈ 0
    end
end

@testset "faulty circuits" begin
    @testset "indeterminate output" begin
        circ = @circuit begin
            r = resistor(0)
            probe = currentprobe(), [+] ⟷ r[1], [-] ⟷ r[2]
        end
        @static if VERSION ≥ v"0.7.0-DEV.2988"
            Test.@test_logs (:warn, "Model output depends on indeterminate quantity") DiscreteModel(circ, 1)
        else
            @test_warn "output depends on indeterminate quantity" DiscreteModel(circ, 1)
        end
    end

    @testset "no solution" begin
        circ = @circuit begin
            d = diode()
            src = currentsource(), [+] ↔ d[+], [-] ↔ d[-]
            probe = voltageprobe(), [+] == d[+], [-] == d[-]
        end
        model = DiscreteModel(circ, 1)
        @test ACME.nn(model) == 1
        y = run!(model, [1.0 1.0])
        @test size(y) == (1, 2)
        @test y[1,1] == y[1,2]
        @test_throws ErrorException run!(model, hcat([Inf]))
        @static if VERSION ≥ v"0.7.0-DEV.2988"
            @test(size(Test.@test_logs((:warn, "Failed to converge while solving non-linear equation."), run!(model, hcat([-1.0])))) == (1, 1))
        else
            @test_warn("Failed to converge", @test size(run!(model, hcat([-1.0]))) == (1, 1))
        end
    end
end

@testset "KDTree" begin
    @testset "size = 4×$num" for num = 1:50
        let ps = rand(4, num)
            t = ACME.KDTree(ps)
            for i in 1:size(ps)[2]
                idx = ACME.indnearest(t, ps[:,i])
                @test ps[:,i] == ps[:,idx]
            end
        end
    end

    @testset "size = 5×10000" begin
        ps = rand(6, 10000)
        t = ACME.KDTree(ps)
        p = rand(6)
        if VERSION ≥ v"0.7.0-DEV.4064"
            best_p = ps[:,argmin(vec(sum(abs2, ps .- p, dims=1)))]
        else
            best_p = ps[:,argmin(vec(sum(abs2, ps .- p, 1)))]
        end
        idx = ACME.indnearest(t, p)
        @test sum(abs2, p - best_p) ≈ sum(abs2, p - ps[:, idx])
    end
end

@testset "HomotopySolver" begin
    nleq = ACME.ParametricNonLinEq((res, J, scratch, z) ->
    let p=scratch[1], Jp=scratch[2]
        res[1] = z[1]^2 - 1 + p[1]
        J[1,1] = 2*z[1]
        Jp[1,1] = 1
    end, 1, 1)
    solver = ACME.HomotopySolver{ACME.SimpleSolver}(nleq, [0.0], [1.0])
    ACME.solve(solver, [-0.5 + rand()])
    @test ACME.hasconverged(solver)
    ACME.solve(solver, [1.5 + rand()])
    @test !ACME.hasconverged(solver)
end

@testset "gensolve/rank_factorize" begin
    a = Rational{BigInt}[1 1 1; 1 1 2; 1 2 1; 1 2 2; 2 1 1; 2 1 2]
    b = Rational{BigInt}[1 2 3 4 5 6; 6 5 4 3 2 1; 1 0 1 0 1 0]
    nullspace = ACME.gensolve(sparse(a'), spzeros(Rational{BigInt}, size(a, 2), 0))[2]
    @test nullspace'*a == spzeros(3, 3)
    c, f = ACME.rank_factorize(sparse(a * b))
    @test c*f == a*b
end

@testset "redcue_pdims" begin
    a = Rational{BigInt}[-1 -1 -4 -3 0 -1; 2 -1 -5 3 -4 0; -2 2 -5 -2 5 1;
                         -5 4 -3 0 5 -5; 4 3 0 -1 0 2; 0 -3 -4 -4 -3 4]
    b = hcat(Rational{BigInt}[1; 2; 3; -2; -1; 0])
    c = Rational{BigInt}[4 2 -1; -1 -3 0; -3 5 3; 0 0 0; -4 -1 -1; -1 -1 5]
    dy = Rational{BigInt}[1 2 3 -2 -1 0]
    ey = hcat(Rational{BigInt}[5])
    fy = Rational{BigInt}[-2 -1 3]
    p = Rational{BigInt}[1 1 1; 1 1 2; 1 2 1; 1 2 2; 2 1 1; 2 1 2]
    dq = Rational{BigInt}[1 2 3 4 5 6; 6 5 4 3 2 1; 1 0 1 0 1 0]
    eq = hcat(Rational{BigInt}[1; 2; 3])
    fq = Rational{BigInt}[1 0 0; 10 0 0; 0 1 0; 0 10 0; 0 0 1; 0 0 10]

    for zxin in (zeros(Rational{BigInt}, 3, 6), Rational{BigInt}[1 2 0 0 2 1; 0 1 2 2 0 1; 0 0 1 0 1 1]),
        zuin in (zeros(Rational{BigInt}, 3, 1), hcat(Rational{BigInt}[1; 2; -1]))
        dq_full = p * dq + fq * zxin
        eq_full = p * eq + fq * zuin
        mats = Dict{Symbol,Array}(:a => a, :b => b, :c => c, :dy => dy, :ey => ey,
            :fy => fy, :dq_full => dq_full, :eq_full => eq_full, :fq => fq)
        mats[:dq_fulls]=Matrix[mats[:dq_full]]
        mats[:eq_fulls]=Matrix[mats[:eq_full]]
        mats[:fqprev_fulls]=Matrix[mats[:eq_full]]
        mats[:fqs]=Matrix[mats[:fq]]
        ACME.reduce_pdims!(mats)
        @test size(mats[:pexps][1], 2) == 3
        @test mats[:pexps][1] * mats[:dqs][1] == mats[:dq_fulls][1]
        @test mats[:pexps][1] * mats[:eqs][1] == mats[:eq_fulls][1]
        zx = (fq'fq) \ fq' * (dq_full - mats[:dq_fulls][1])
        zu = (fq'fq) \ fq' * (eq_full - mats[:eq_fulls][1])
        @test mats[:a] == a - c*zx
        @test mats[:b] == b - c*zu
        @test mats[:dy] == dy - fy*zx
        @test mats[:ey] == ey - fy*zu
    end

    # TODO: Also check dq, eq, fqprev for a test case with desomposed non-linearity
end

@testset "nonlinearity decomposition" begin
    circ = @circuit begin
        src1 = voltagesource()
        probe1 = currentprobe()
        d1 = diode(), [+] ⟷ src1[+]
        d2 = diode(), [+] ⟷ d1[-], [-] ⟷ probe1[+]
        probe1[-] ⟷ src1[-]
    end
    add!(circ, :src2, voltagesource())
    add!(circ, :probe2, currentprobe())
    add!(circ, :d3, diode())
    connect!(circ, (:src2, :+), (:d3, :+))
    connect!(circ, (:d3, :-), (:probe2, :+))
    connect!(circ, (:probe2, :-), (:src2, :-))
    model = DiscreteModel(circ, 1, decompose_nonlinearity=false)
    y = run!(model, hcat([2.0; 1.0]))
    @test ACME.nn(model, 1) == 3
    @test y[1] ≈ 1e-12*(exp(1/25e-3)-1)
    @test y[2] ≈ 1e-12*(exp(1/25e-3)-1)
    model = DiscreteModel(circ, 1)
    y = run!(model, hcat([2.0; 1.0]))
    # single diode is extracted first, although it was added last
    @test ACME.nn(model, 1) == 1
    @test ACME.nn(model, 2) == 2
    @test y[1] ≈ y[2] ≈ 1e-12*(exp(1/25e-3)-1)
end

@testset "sources/probes" begin
# sources and probes with internal resistance/conductance
    circ = @circuit begin
        src = currentsource(100e-3, gp=1//100000)
        probe = voltageprobe(), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, zeros(0,1)) ≈ [100000*100e-3]

    circ = @circuit begin
        src = currentsource(gp=1//100000)
        probe = voltageprobe(), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, hcat([100e-3])) ≈ [100000*100e-3]

    circ = @circuit begin
        src = currentsource(100e-3)
        probe = voltageprobe(gp=1//100000), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, zeros(0,1)) ≈ [100000*100e-3]

    circ = @circuit begin
        src = voltagesource(10, rs=100000)
        probe = currentprobe(), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, zeros(0,1)) ≈ [10/100000]

    circ = @circuit begin
        src = voltagesource(rs=100000)
        probe = currentprobe(), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, hcat([10.0])) ≈ [10/100000]

    circ = @circuit begin
        src = voltagesource(10)
        probe = currentprobe(rs=100000), [+] ⟷ src[+], [-] ⟷ src[-]
    end
    model = DiscreteModel(circ, 1)
    @test run!(model, zeros(0,1)) ≈ [10/100000]
end

@testset "BJT" begin
    isc=1e-6
    ise=2e-6
    ηc=1.1
    ηe=1.0
    βf=100
    βr=10
    @testset "Ebers-Moll $typ" for (typ, ib) in ((:npn, 1e-3), (:pnp, -1e-3))
        circ = @circuit begin
            t = bjt(typ, isc=isc, ise=ise, ηc=ηc, ηe=ηe, βf=βf, βr=βr)
            isrc = currentsource(), [+] ⟷ t[base]
            vsrc = voltagesource(), [-] ⟷ isrc[-]
            veprobe = voltageprobe(), [+] ⟷ t[base], [-] ⟷ isrc[-]
            vcprobe = voltageprobe(), [+] ⟷ t[base], [-] ⟷ vsrc[+]
            ieprobe = currentprobe(), [+] ⟷ t[emitter], [-] ⟷ isrc[-]
            icprobe = currentprobe(), [+] ⟷ t[collector], [-] ⟷ vsrc[+]
        end
        model = DiscreteModel(circ, 1)
        N = 100
        output = run!(model, [range(0, stop=ib, length=N)'; range(1, stop=-1, length=N÷2)' range(-1, stop=1, length=N÷2)'])
        if typ == :pnp
            output = -output
        end
        for n in 1:N
            ve, vc, ie, ic = output[:,n]
            @test isapprox(ie, ise*(exp(ve/(ηe*25e-3))-1) - βr/(1+βr)*isc*(exp(vc/(ηc*25e-3))-1), atol=1e-10)
            @test isapprox(ic, -βf/(1+βf)*ise*(exp(ve/(ηe*25e-3))-1) + isc*(exp(vc/(ηc*25e-3))-1), atol=1e-10)
        end
    end
    ηcl=1.2
    ηel=1.3
    prog = Progress(2^9, "Testing Gummel-Poon model: ")
    @testset "Gummel-Poon $typ" for ile in (0, 50e-9), ilc in (0, 100e-9),
            ηcl in (ηc, 1.2), ηel in (ηe, 1.1),
            vaf in (Inf, 10), var in (Inf, 50),
            ikf in (Inf, 50e-3), ikr in (Inf, 500e-3),
            (typ, ib) in ((:npn, 1e-3), (:pnp, -1e-3))
        circ = @circuit begin
            t = bjt(typ, isc=isc, ise=ise, ηc=ηc, ηe=ηe, βf=βf, βr=βr, ile=ile,
                    ilc=ilc, ηcl=ηcl, ηel=ηel, vaf=vaf, var=var, ikf=ikf, ikr=ikr)
            isrc = currentsource(), [+] ⟷ t[base]
            vsrc = voltagesource(), [-] ⟷ isrc[-]
            veprobe = voltageprobe(), [+] ⟷ t[base], [-] ⟷ isrc[-]
            vcprobe = voltageprobe(), [+] ⟷ t[base], [-] ⟷ vsrc[+]
            ieprobe = currentprobe(), [+] ⟷ t[emitter], [-] ⟷ isrc[-]
            icprobe = currentprobe(), [+] ⟷ t[collector], [-] ⟷ vsrc[+]
        end
        model = DiscreteModel(circ, 1)
        N = 100
        output = run!(model, [range(0, stop=ib, length=N)'; range(1, stop=-1, length=N÷2)' range(-1, stop=1, length=N÷2)'])
        if typ == :pnp
            output = -output
        end
        for n in 1:N
            ve, vc, ie, ic = output[:,n]
            i_f = βf/(1+βf)*ise*(exp(ve/(ηe*25e-3))-1)
            i_r = βr/(1+βr)*isc*(exp(vc/(ηc*25e-3))-1)
            icc = (2*(1-ve/var-vc/vaf))/(1+sqrt(1+4(i_f/ikf+i_r/ikr))) * (i_f - i_r)
            ibe = 1/βf*i_f + ile*(exp(ve/(ηel*25e-3))-1)
            ibc = 1/βr*i_r + ilc*(exp(vc/(ηcl*25e-3))-1)
            @test isapprox(ie, icc + ibe, atol=1e-10)
            @test isapprox(ic, -icc + ibc, atol=1e-10)
        end
        next!(prog)
    end

    rb=100
    re=10
    rc=20
    @testset "internal resistances ($typ)" for (typ, ib, vce) in ((:npn, 1e-3, 1), (:pnp, -1e-3, -1))
        circ = @circuit begin
            t1 = bjt(typ)
            rbref = resistor(rb)
            rcref = resistor(rc)
            reref = resistor(re)
            isrc1 = currentsource(ib)
            vscr1 = voltagesource(vce)
            veprobe1 = voltageprobe()
            vcprobe1 = voltageprobe()
            ieprobe1 = currentprobe()
            icprobe1 = currentprobe()
            t1[base] ⟷ rbref[1]
            rbref[2] ⟷ isrc1[+] ⟷ veprobe1[+] ⟷ vcprobe1[+]
            t1[collector] ⟷ rcref[1]
            rcref[2] ⟷ icprobe1[+]
            vcprobe1[-] ⟷ icprobe1[-] ⟷ vscr1[+]
            t1[emitter] ⟷ reref[1]
            reref[2] ⟷ ieprobe1[+]
            veprobe1[-] ⟷ ieprobe1[-] ⟷ vscr1[-] ⟷ isrc1[-]
            t2 = bjt(typ, rb=rb, re=re, rc=rc)
            isrc2 = currentsource(ib)
            vscr2 = voltagesource(vce)
            veprobe2 = voltageprobe()
            vcprobe2 = voltageprobe()
            ieprobe2 = currentprobe()
            icprobe2 = currentprobe()
            t2[base] ⟷ isrc2[+] ⟷ veprobe2[+] ⟷ vcprobe2[+]
            t2[collector] ⟷ icprobe2[+]
            vcprobe2[-] ⟷ icprobe2[-] ⟷ vscr2[+]
            t2[emitter] ⟷ ieprobe2[+]
            veprobe2[-] ⟷ ieprobe2[-] ⟷ vscr2[-] ⟷ isrc2[-]
        end
        model = DiscreteModel(circ, 1)
        output = run!(model, zeros(0,1))
        @test output[1:4,:] ≈ output[5:8,:]
    end
end

function checksteady!(model)
    x_steady = steadystate!(model)
    for s in model.solvers
        ACME.set_resabstol!(s, 1e-13)
    end
    run!(model, zeros(1, 1))
    return model.x ≈ x_steady
end

function linearization_error!(model, amplitude)
    linmodel = linearize(model)
    N = 50000
    u = [amplitude * sin(π/2 * n^2/N) for n in 0:N]'
    steadystate!(model)
    steadystate!(linmodel)
    y = run!(model, u)
    ylin = run!(linmodel, u)
    return maximum(abs, y-ylin)
end

@testset "examples" begin
    @testset "sallenkey" begin
        include(joinpath(dirname(@__FILE__), "..", "examples", "sallenkey.jl"))
         model=sallenkey()
         println("Running sallenkey")
         y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
         @test size(y) == (1,44100)
         # TODO: further validate y
         @test checksteady!(model)
     end

     @testset "diodeclipper" begin
        include(joinpath(dirname(@__FILE__), "..", "examples", "diodeclipper.jl"))
        model=diodeclipper()
        println("Running diodeclipper")
        @test ACME.np(model, 1) == 1
        y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
        @test checksteady!(model)

        @test linearization_error!(model, 1e-3) < 1e-15

        model = diodeclipper(solver=HomotopySolver{SimpleSolver})
        runner = ModelRunner(model, false)
        u = sin.(2π*1000/44100*(0:44099)')
        y = run!(runner, u)
        run!(runner, y, u)
        alloc = @allocated run!(runner, y, u)
        if alloc > 0
            warn("Allocated $alloc in run! of HomotopySolver{SimpleSolver} base ModelRunner")
        end
    end

    @testset "birdie" begin
        include(joinpath(dirname(@__FILE__), "..", "examples", "birdie.jl"))
        model=birdie(vol=0.8)
        ACME.solve(model.solvers[1], [0.003, -0.0002])
        @assert all(ACME.hasconverged, model.solvers)
        println("Running birdie with fixed vol")
        @test ACME.np(model, 1) == 2
        y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
        @test checksteady!(model)

        @test linearization_error!(model, 1e-4) < 1e-7

        model=birdie()
        println("Running birdie with varying vol")
        @test ACME.np(model, 1) == 3
        y = run!(model, [sin.(2π*1000/44100*(0:44099)'); range(1, stop=0, length=44100)']; showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
    end

    @testset "superover" begin
        include(joinpath(dirname(@__FILE__), "..", "examples", "superover.jl"))
        model=superover(drive=1.0, tone=1.0, level=1.0)
        println("Running superover with fixed potentiometer values")
        @test ACME.np(model, 1) == 5
        y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
        @test checksteady!(model)
        @test linearization_error!(model, 1e-4) < 1e-4 # SuperOver really is not very linear...

        circ=superover(Circuit, drive=1.0, tone=1.0, level=1.0)
        println("Running simplified superover with fixed potentiometer values")
        add!(circ, :vbsrc, voltagesource(4.5))
        connect!(circ, (:vbsrc, :+), :vb)
        connect!(circ, (:vbsrc, :-), :gnd)
        model = DiscreteModel(circ, 1/44100)
        @test ACME.np(model, 1) == 2
        @test ACME.np(model, 2) == 1
        @test ACME.np(model, 3) == 2
        y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
        @test_broken checksteady!(model)
        @test_broken linearization_error!(model, 1e-4) < 1e-4 # SuperOver really is not very linear...

        println("Running simplified, non-decomposed superover with fixed potentiometer values")
        model = DiscreteModel(circ, 1/44100, decompose_nonlinearity=false)
        @test ACME.np(model, 1) == 5
        y = run!(model, sin.(2π*1000/44100*(0:44099)'); showprogress=false)
        @test size(y) == (1,44100)
        # TODO: further validate y
        @test checksteady!(model)
        @test linearization_error!(model, 1e-4) < 1e-4 # SuperOver really is not very linear...

        model=superover()
        println("Running superover with varying potentiometer values")
        @test ACME.np(model, 1) == 11
        y = run!(model, [sin.(2π*1000/44100*(0:999)'); range(1, stop=0, length=1000)'; range(0, stop=1, length=1000)'; range(1, stop=0, length=1000)']; showprogress=false)
        @test size(y) == (1,1000)
        # TODO: further validate y

        circ=superover(Circuit)
        println("Running simplified superover with varying potentiometer values")
        add!(circ, :vbsrc, voltagesource(4.5))
        connect!(circ, (:vbsrc, :+), :vb)
        connect!(circ, (:vbsrc, :-), :gnd)
        model = DiscreteModel(circ, 1/44100)
        @test ACME.np(model, 1) == 2
        @test ACME.np(model, 2) == 2
        @test ACME.np(model, 3) == 2
        @test ACME.np(model, 4) == 4
        y = run!(model, [sin.(2π*1000/44100*(0:999)'); range(1, stop=0, length=1000)'; range(0, stop=1, length=1000)'; range(1, stop=0, length=1000)']; showprogress=false)
        @test size(y) == (1,1000)
        # TODO: further validate y
    end
end


text:
struct LRFinderResult{T}
    losses::Vector{T}
    lrs::Vector{T}
    steepest::T
    mindiv10::T
end

function LRFinderResult(losses::Vector{T}, lrs) where T
    return LRFinderResult(losses, lrs, zero(T), zero(T))
end

mutable struct LRFinderPhase <: FluxTraining.Phases.AbstractTrainingPhase
    startlr
    endlr
    steps
    β
    divergefactor
    result
end

function LRFinderPhase(; startlr = 1e-7, endlr = 10, steps = 100, β = 0.02, divergefactor = 4, result = nothing)
    return LRFinderPhase(startlr, endlr, steps, β, divergefactor, result)
end

function FluxTraining.fitepochphase!(
        learner::Learner,
        phase::LRFinderPhase)

    dataiter = FluxTraining.getdataiter(TrainingPhase(), learner)
    if dataiter === nothing
        throw(CancelEpochException("No data found for phase $(typeof(phase))"))
    end

    metric = SmoothLoss(phase.β)
    metricscb = Metrics(metric)
    losses = Float32[]
    lrs = Float32[]

    withfields(
        learner,
        model = (FluxTraining.model!, deepcopy(cpu(learner.model))),
        optimizer = deepcopy(learner.optimizer)
        ) do

        withcallbacks(learner, metricscb, Scheduler(Dict())) do
            #schedule = Schedule([0, phase.steps], [phase.startlr, phase.endlr], Animations.expin(100000000))
            #setschedules!(learner, phase, LearningRate => schedule)
            bestloss = Inf

            handle(EpochBegin(), learner, phase)
            for (i, batch) in zip(1:phase.steps, dataiter)
                lr = phase.startlr * (phase.endlr / phase.startlr) ^ (i / phase.steps)
                learner.optimizer.eta = lr
                FluxTraining.fitbatch!(learner, batch, phase)

                loss = FluxTraining.stepvalue(metric)
                bestloss = min(bestloss, loss)
                push!(losses, loss)
                push!(lrs, lr)

                # Stop if loss is diverging
                if loss > bestloss * phase.divergefactor
                    break
                end
            end
            handle(EpochEnd(), learner, phase)
        end
    end

    phase.result = LRFinderResult(losses, lrs)
end

function plotlrfind(lrs, losses)
    f = Figure(resolution = (700, 400))
    f[1, 1] = ax = Axis(f)

    xtickvalues = (10.) .^(-10:10)
    ax.xlabel = "Learning rate"
    ax.xticks = log.(xtickvalues)
    ax.xtickformat[] = vals -> string.(xtickvalues)
    ax.xticklabelsize = 10.

    ax.ylabel = "Loss"
    ax.yticks = []
    ax.ytickformat[] = vals -> string.(log.(vals))
    ax.yticklabelsize = 10.

    AbstractPlotting.lines!(ax, log.(lrs), log.(losses), axis = (xticks = LinearTicks(10),))
    return f
end

plotlrfind(phase::LRFinderPhase) = plotlrfind(phase.result)
plotlrfind(lrresult::LRFinderResult) = plotlrfind(lrresult.lrs, lrresult.losses)


text:
# Atomic Functions

# TODO: support for 64-bit atomics via atom_cmpxchg (from cl_khr_int64_base_atomics)

# "atomic operations on 32-bit signed, unsigned integers and single precision
#  floating-point to locations in __global or __local memory"

const atomic_integer_types = [UInt32, Int32]
# TODO: 64-bit atomics with ZE_DEVICE_MODULE_FLAG_INT64_ATOMICS
# TODO: additional floating-point atomics with ZE_extension_float_atomics
const atomic_memory_types = [AS.Local, AS.Global]

# generically typed

for gentype in atomic_integer_types, as in atomic_memory_types
@eval begin

@device_function atomic_add!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_add", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_sub!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_sub", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_inc!(p::LLVMPtr{$gentype,$as}) =
    @builtin_ccall("atomic_inc", $gentype, (LLVMPtr{$gentype,$as},), p)

@device_function atomic_dec!(p::LLVMPtr{$gentype,$as}) =
    @builtin_ccall("atomic_dec", $gentype, (LLVMPtr{$gentype,$as},), p)

@device_function atomic_min!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_min", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_max!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_max", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_and!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_and", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_or!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_or", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_xor!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_xor", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_xchg!(p::LLVMPtr{$gentype,$as}, val::$gentype) =
    @builtin_ccall("atomic_xchg", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype), p, val)

@device_function atomic_cmpxchg!(p::LLVMPtr{$gentype,$as}, cmp::$gentype, val::$gentype) =
    @builtin_ccall("atomic_cmpxchg", $gentype,
                   (LLVMPtr{$gentype,$as}, $gentype, $gentype), p, cmp, val)

end
end

# specifically typed

for as in atomic_memory_types
@eval begin

@device_function atomic_xchg!(p::LLVMPtr{Float32,$as}, val::Float32) =
    @builtin_ccall("atomic_xchg", Float32, (LLVMPtr{Float32,$as}, Float32,), p, val)

# XXX: why is only xchg supported on floats? isn't it safe for cmpxchg too,
#      which should only perform bitwise comparisons?
@device_function atomic_cmpxchg!(p::LLVMPtr{Float32,$as}, cmp::Float32, val::Float32) =
    reinterpret(Float32, atomic_cmpxchg!(reinterpret(LLVMPtr{UInt32,$as}, p),
                                         reinterpret(UInt32, cmp),
                                         reinterpret(UInt32, val)))

end
end

# documentation

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `old + val` and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_add!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `old - val` and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_sub!

"""
Swaps the old value stored at location `p` with new value given by `val`.
Returns old value.
"""
atomic_xchg!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute (`old` + 1) and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_inc!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute (`old` - 1) and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_dec!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `(old == cmp) ? val : old` and store result at location pointed by `p`.
The function returns `old`.
"""
atomic_cmpxchg!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `min(old, val)` and store minimum value at location pointed by `p`. The
function returns `old`.
"""
atomic_min!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `max(old, val)` and store maximum value at location pointed by `p`. The
function returns `old`.
"""
atomic_max

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `old & val` and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_and!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `old | val` and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_or!

"""
Read the 32-bit value (referred to as `old`) stored at location pointed by `p`.
Compute `old ^ val` and store result at location pointed by `p`. The function
returns `old`.
"""
atomic_xor!

#
# High-level interface
#

# prototype of a high-level interface for performing atomic operations on arrays
#
# this design could be generalized by having atomic {field,array}{set,ref} accessors, as
# well as acquire/release operations to implement the fallback functionality where any
# operation can be applied atomically.

const inplace_ops = Dict(
    :(+=) => :(+),
    :(-=) => :(-),
    :(*=) => :(*),
    :(/=) => :(/),
    :(&=) => :(&),
    :(|=) => :(|),
    :(⊻=) => :(⊻)
)

struct AtomicError <: Exception
    msg::AbstractString
end

Base.showerror(io::IO, err::AtomicError) =
    print(io, "AtomicError: ", err.msg)

"""
    @atomic a[I] = op(a[I], val)
    @atomic a[I] ...= val

Atomically perform a sequence of operations that loads an array element `a[I]`, performs the
operation `op` on that value and a second value `val`, and writes the result back to the
array. This sequence can be written out as a regular assignment, in which case the same
array element should be used in the left and right hand side of the assignment, or as an
in-place application of a known operator. In both cases, the array reference should be pure
and not induce any side-effects.

!!! warn
    This interface is experimental, and might change without warning.  Use the lower-level
    `atomic_...!` functions for a stable API.
"""
macro atomic(ex)
    # decode assignment and call
    if ex.head == :(=)
        ref = ex.args[1]
        rhs = ex.args[2]
        Meta.isexpr(rhs, :call) || throw(AtomicError("right-hand side of an @atomic assignment should be a call"))
        op = rhs.args[1]
        if rhs.args[2] != ref
            throw(AtomicError("right-hand side of a non-inplace @atomic assignment should reference the left-hand side"))
        end
        val = rhs.args[3]
    elseif haskey(inplace_ops, ex.head)
        op = inplace_ops[ex.head]
        ref = ex.args[1]
        val = ex.args[2]
    else
        throw(AtomicError("unknown @atomic expression"))
    end

    # decode array expression
    Meta.isexpr(ref, :ref) || throw(AtomicError("@atomic should be applied to an array reference expression"))
    array = ref.args[1]
    indices = Expr(:tuple, ref.args[2:end]...)

    esc(quote
        $atomic_arrayset($array, $indices, $op, $val)
    end)
end

# FIXME: make this respect the indexing style
@inline atomic_arrayset(A::AbstractArray{T}, Is::Tuple, op::Function, val) where {T} =
    atomic_arrayset(A, Base._to_linear_index(A, Is...), op, convert(T, val))

# native atomics
for (op,impl) in [(+)      => atomic_add!,
                  (-)      => atomic_sub!,
                  (&)      => atomic_and!,
                  (|)      => atomic_or!,
                  (⊻)      => atomic_xor!,
                  Base.max => atomic_max!,
                  Base.min => atomic_min!]
    @eval @inline atomic_arrayset(A::oneDeviceArray{T}, I::Integer, ::typeof($op),
                                  val::T) where {T <: Union{Int32,UInt32}} =
        $impl(pointer(A, I), val)
end

# fallback using compare-and-swap
function atomic_arrayset(A::AbstractArray{T}, I::Integer, op::Function, val) where {T}
    ptr = pointer(A, I)
    old = Base.unsafe_load(ptr, 1)
    while true
        cmp = old
        new = convert(T, op(old, val))
        old = atomic_cmpxchg!(ptr, cmp, new)
        (old == cmp) && return new
    end
end


text:
#Demo 5: Heat map
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Example using ColorSchemes
#(run `using ColorSchemes` before running demo to test)
if @isdefined(ColorSchemes)
cs = ColorSchemes.jet1
colorscale = InspectDR.ColorScale(get(cs, range(0, stop=1, length=100)))
end

#Input data
#-------------------------------------------------------------------------------
x = collect(0:1.0:10.0)
y = collect(0:1.0:4.0)
zi = float((1:10) * reshape(1:4, 1, :))
zmax = maximum(zi); zmin = minimum(zi)
z = (zi.+.25*zmax).*1e9
lz = (zi.-zmin.+1e-3)
lz = log10.(lz.*lz.*1e3)

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Heat Maps")
mplot.layout[:ncolumns] = 1

plot = add(mplot, InspectDR.Plot2D(:lin, [:lin, :lin],
		title="Heatmaps (upper: linear, lower: log)", xlabel="x", ylabels=["y", "y"])
	)
	plot.layout[:enable_colorscale] = true
	strip_linz, strip_logz = plot.strips

if @isdefined(ColorSchemes)
	strip_linz.colorscale = colorscale
	strip_logz.colorscale = colorscale
end

	#Define extents & scales:
#	strip_logz.grid = InspectDR.GridRect(vmajor=true, hmajor=true, hminor=true)
	wfrm = addheatmap(plot, x, y, z, strip=1)
	wfrm = addheatmap(plot, x, y, lz, strip=2)

gplot = display(InspectDR.GtkDisplay(), mplot)

:DONE


text:
abstract type LaneFollowingDriver <: DriverModel{LaneFollowingAccel} end
track_longitudinal!(model::LaneFollowingDriver, v_ego::Float64, v_oth::Float64, headway::Float64) = model # do nothing by default

function observe!(model::LaneFollowingDriver, scene::Frame{Entity{S, D, I}}, roadway::Roadway, egoid::I) where {S, D, I}

    vehicle_index = findfirst(egoid, scene)

    fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())

    v_ego = vel(scene[vehicle_index].state)
    v_oth = NaN
    headway = NaN

    if fore.ind != nothing
        v_oth = vel(scene[fore.ind].state)
        headway = fore.Δs
    end

    track_longitudinal!(model, v_ego, v_oth, headway)

    return model
end

mutable struct StaticLaneFollowingDriver <: LaneFollowingDriver
    a::LaneFollowingAccel
end
StaticLaneFollowingDriver() = StaticLaneFollowingDriver(LaneFollowingAccel(0.0))
StaticLaneFollowingDriver(a::Float64) = StaticLaneFollowingDriver(LaneFollowingAccel(a))
get_name(::StaticLaneFollowingDriver) = "ProportionalSpeedTracker"
Base.rand(rng::AbstractRNG, model::StaticLaneFollowingDriver) = model.a
Distributions.pdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : 0.0
Distributions.logpdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : -Inf


text:
using Polyhedra
using CDDLib
using MAT

include("../../src/Reachability/ReachabilityPolyhedra.jl")
using ReachabilityPolyhedra

function compute(algorithm="reach_homog")

# system's matrix
A = [0. 1. 0. 0. 0. 0. ;
-0.417533 -3.1931759963 39.24 0. -14.825331 11.123344 ;
0. 0. 0. 1. 0. 0. ;
0.0417533 0.31931759963 -4.905 0. 1.4825331 -1.1123344 ;
0.0638407957 -0.32473339016573 0. 0. -3.7332068901 -0.7007592976 ;
0.0853437452 -0.72366802635628 0. 0. -5.9714023436 -2.2736115136]

# number of variables (dimension)
n = size(A, 1)

# initial set of states
file = matopen("X0.mat")

#=
F = read(file, "F");
g = read(file, "g");
g = g[:];  # cast into a column vector
# this is a half-space representation. use the methods X0.A and X0.b to access
# the matrix/vectors such that X0 = {x : Ax <= b}
X0HRep = SimpleHRepresentation(F, g);
X0 = polyhedron(X0HRep, CDDLibrary());
=#

# initial states : hyperrectangle
g1 = read(file, "g1")[:]
g2 = read(file, "g1")[:]
g3 = read(file, "g1")[:]
F1 = read(file, "F1")
F2 = read(file, "F2")
F3 = read(file, "F3")

X0hat = Array{CDDPolyhedron}(3)
X0hat[1] = polyhedron(SimpleHRepresentation(F1, g1))
X0hat[2] = polyhedron(SimpleHRepresentation(F2, g2))
X0hat[3] =  polyhedron(SimpleHRepresentation(F3, g3))

# number of samples
N = 250

# ==== case with inputs ===

# this polyhedron is "null" in n dims
#nullPoly_6d = polyhedron(SimpleVRepresentation([0.0 0.0 0.0 0.0 0.0 0.0]));
# we specify the backend using a second argument, as CDDLibrary()
# otherwise, the object will be of type Polyhedra.SimplePolyhedron{6,Float64},
# and we cannot do operations with it
nullPoly_6d = polyhedron(SimpleVRepresentation(zeros(1, 6)), CDDLibrary());
U = [nullPoly_6d for i in 1:N]

# === solve ===

if algorithm == "reach_homog"
    RSet = reach_homog(A, X0hat, N)
elseif algorithm == "reach"
    RSet = reach(A, X0hat, N, U)
else
    error("algorithm unkwnown")
end

end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
__precompile__()

module DASSL

export dasslIterator, dasslSolve

using Reexport
@reexport using DiffEqBase
using LinearAlgebra

import DiffEqBase: solve

export dassl

include("common.jl")

const MAXORDER = 6
const MAXIT = 10

mutable struct JacData
    a::Real
    jac # Jacobian matrix for the newton solver
end

function dasslStep(channel,
                   F,
                   y0::AbstractVector{T},
                   tstart::Real;
                   reltol   = 1e-3,
                   abstol   = 1e-5,
                   initstep = 1e-4,
                   maxstep  = Inf,
                   minstep  = 0,
                   maxorder = MAXORDER,
                   dy0      = zero(y0),
                   tstop    = Inf,
                   norm     = dassl_norm,
                   weights  = dassl_weights,
                   factorize_jacobian = true, # whether to store factorized version of jacobian
                   jacobian = numerical_jacobian(F,reltol,abstol,weights), # computes the quantity F/dy+a*dF/dy'
                   args...) where T<:Number

    n  = length(y0)
    jd = JacData(zero(tstart),zeros(T,n,n)) # generate a dummy
                                   # jacobian, it will be replaced by
                                   # a real one after running stepper

    ord      = 1                # initial method order
    tout     = [tstart]         # initial time
    h        = initstep         # current step size
    yout     = Array{typeof(y0)}(undef,1)
    yout[1]  = y0
    dyout    = Array{typeof(y0)}(undef,1)
    dyout[1] = dy0              # initial guess for dy0
    num_rejected = 0            # number of rejected steps
    num_fail = 0                # number of consecutive failures

    # wrapper around the jacobian function
    if factorize_jacobian
        computejac = (t,y,dy,a) -> LinearAlgebra.factorize(jacobian(t,y,dy,a))
    else
        computejac = jacobian
    end

    # This is the trick to improve on the initial guess for
    # y0.  Basically we run a one iteration of stepper and use the
    # result as the initial derivative
    (_,_,_,dyout[1],_)=stepper(1,tout,yout,dyout,10*eps(one(tstart)),F,jd,computejac,weights(y0,1,1),v->norm(v,weights(y0,1,1)),1)

    while tout[end] < tstop

        hmin = max(4*eps(one(tstart)),minstep)
        h = min(h,maxstep,tstop-tout[end])

        if h < hmin
            info("Stepsize too small (h=$h at t=$(tout[end]).")
            break
        elseif num_fail >= -2/3*log(eps(one(tstart)))
            info("Too many ($num_fail) failed steps in a row (h=$h at t=$(tout[end]).")
            break
        end

        # error weights
        wt = weights(yout[end],reltol,abstol)
        normy(v) = norm(v,wt)

        (status,err,yn,dyn,jd)=stepper(ord,tout,yout,dyout,h,F,jd,computejac,wt,normy,maxorder)

        if status < 0
            # Early failure: Newton iteration failed to converge, reduce
            # the step size and try again

            # increase the failure counter
            num_fail     += 1
            # keep track of the total number of rejected steps
            num_rejected += 1
            # reduce the step by 25%
            h *= 1/4
            continue

        elseif err > 1
            # local error is too large.  Step is rejected, and we try
            # again with new step size and order.

            # increase the failure counter
            num_fail     += 1
            # keep track of the total number of rejected steps
            num_rejected += 1
            # temporarily push the new step to the yout (needed by newStepOrder)
            push!(tout,tout[end]+h); push!(yout,yn)
            # determine the new step size and order, excluding the current step
            (r,ord) = newStepOrder(tout,yout,normy,err,ord,num_fail,maxorder)

            # pop the temporary steps
            pop!(tout); pop!(yout)
            h *= r
            continue

        else
            ####################
            # step is accepted #
            ####################

            # reset the failure counter
            num_fail      = 0

            # save the results
            push!(tout, tout[end]+h)
            push!(yout, yn)
            push!(dyout,dyn)

            # remove old results
            if length(tout) > ord+3
                popfirst!(tout)
                popfirst!(yout)
                popfirst!(dyout)
            end

            push!(channel,(tout[end],yout[end],dyout[end]))

            # determine the new step size and order, including the current step
            (r,ord) = newStepOrder(tout,yout,normy,err,ord,num_fail,maxorder)

            h *= r
        end

    end

end

# the iterator version of the dasslSolve
dasslIterator(F, y0, t0; args...) = Channel((channel)->dasslStep(channel,F, y0, t0; args...))

# solves the equation F with initial data y0 over for times t in tspan=[t0,t1]
function dasslSolve(F, y0::AbstractVector, tspan; dy0 = zero(y0), args...)
    tout  = Array{typeof(tspan[1])}(undef,1)
    yout  = Array{typeof(y0)}(undef,1)
    dyout = Array{typeof(y0)}(undef,1)
    tout[1]  = tspan[1]
    yout[1]  = y0
    dyout[1] = dy0
    for (t, y, dy) in dasslIterator(F, y0, tspan[1]; dy0=dy0, tstop=tspan[end], args...)
        push!( tout,  t)
        push!( yout,  y)
        push!(dyout, dy)
        if t >= tspan[end]
            break
        end
    end
    return (tout,yout,dyout)
end

# A scalar version of dasslSolve, implemented as a wrapper around dasslSolve
function dasslSolve(F, y0::Number, tspan; args...)
    (tout,yout,dyout) = dasslSolve(F,[y0],tspan; args...)
    return (tout,map(first,yout),map(first,dyout))
end

function newStepOrder(t::AbstractVector,
                      y::AbstractVector,
                      normy,
                      err,
                      k::Integer,
                      num_fail::Integer,
                      maxorder::Integer)

    if length(t) != length(y)
        error("incompatible size of y and t")
    end

    available_steps = length(t) # including the step t_{n+1}

    if num_fail >= 3
        # probably, the step size was drastically decreased for
        # several steps in a row, so we reduce the order to one and
        # further decrease the step size
        (r,order) = (1/4,1)

    elseif num_fail == 1 && available_steps == 1 && err > 1
        # the last step was accepted, increase order to two
        (r,order) = (1/4,1)

    elseif num_fail == 0 && available_steps == 2 && err < 1
        # The first successfull step, try increasing the order without changing the step size
        (r,order) = (1, 2)

    elseif available_steps < k+2
        # we are at the beginning of the integration, we don't have
        # enough steps to run newStepOrderContinuous, we have to rely
        # on a crude order/stepsize selection
        if num_fail == 0
            # previous step was accepted so we can increase the order
            # and the step size
            # (r,order) = (2.0,min(k+1,maxorder))
            (r,order) = ((2*err+1/10000)^(-1/(k+1)), k)
            # (r,order) = (2.0,1)
        else
            # @todo I am not sure about the choice of order
            #
            # previous step was rejected, we have to decrease the step
            # size and order
            (r,order) = (1/4,1)
        end

    else
        # we have at least k+3 previous steps available, so we can
        # safely estimate the order k-2, k-1, k and possibly k+1
        (r,order) = newStepOrderContinuous(t,y,normy,err,k,maxorder)
        # this function prevents from step size changing too rapidly
        r = normalizeStepSize(r,num_fail)

        # if the previous step failed don't increase the order
        if num_fail > 0
            order = min(order,k)
        end

    end

    return r, order

end

function newStepOrderContinuous(t::AbstractVector,
                                y::AbstractVector,
                                normy,
                                err,
                                k::Integer,
                                maxorder::Integer)

    # compute the error estimates of methods of order k-2, k-1, k and
    # (if possible) k+1
    errors  = errorEstimates(t,y,normy,k)
    errors[k] = err
    ne = length(errors)         # == k or k+1

    # we want to be conservative in our choice of order so we consider
    # not only a monotone sequences but geometrically monotone
    # sequences
    errors_dec  = all(diff([errors[i] for i=max(k-2,1):min(ne,maxorder)]).<0)
    errors_inc  = all(diff([errors[i] for i=max(k-2,1):min(ne,maxorder)]).>0)

    if ne == k+1 && errors_dec
        # If we can estimate the error for k+1 order and the Taylor
        # expansion errors form a decreasing sequence, we can safely
        # increase the order
        order = min(k+1,maxorder)
    elseif ne > 1 && errors_inc
        # Taylor expansion errors form an increasing sequence, we
        # should decrease the order
        order = max(k-1,1)
    else
        # otherwise, leave the current order
        order = k
    end

    # error estimate for the next step
    est = errors[order]

    # initial guess for the new step size multiplier
    r = (2*est+1/10000)^(-1/(order+1))

    return r, order

end

# Based on whether the previous steps were successful we determine
# the new step size
#
# num_fail is the number of steps that failed before this step, r is a
# suggested step size multiplier.
function normalizeStepSize(r, num_fail)

    if num_fail == 0
        # previous step was accepted
        if r >= 2
            r = 2.0
        elseif r < 1
            # choose r from between 0.5 and 0.9
            r = max(1/2,min(r,9/10))
        else
            r = 1.0
        end

    elseif num_fail == 1
        # previous step failed, we slightly decrease the step size,
        # the resulting r is between 0.25 and 0.9
        r = max(1/4,9/10*min(r,1))

    elseif num_fail == 2
        # previous step failed for a second time, error estimates are
        # probably not reliable so decrease the step size
        r = 1/4

    end

    return r

end

# this function estimates the errors of methods of order k-2,k-1,k,k+1
# and returns the estimates as an array.
#
# This method requires at least k previous steps (not including the
# predicted step), so the length of tables t and y should be at least
# k+2. In particular, these estimates won't work for the first step.
#
# The error estimates for order `kest` are computed as the maximal
# derivative `err=norm(h^{(kest+1)}*u^{(kest+1)})` of an interpolating
# polynomial constructed from kest+1 steps at
# t_{n+1},...,t_{n-kest}.  So to estimate the error of e.g. order
# `kest=1` we need steps `t_{n+1},t_{n},t_{n-1}`, in general this
# method needs k+1 previous steps and one future step.
#
# additionally, if the previous k+1 steps were made with the same time
# step, we can estimate the error for order k+1.
#

# here t is an array of times    t=[t_1, ..., t_n, t_{n+1}]
# and y is an array of solutions y=[y_1, ..., y_n, y_{n+1}]
function errorEstimates(t::AbstractVector,
                        y::AbstractVector,
                        normy,
                        k::Integer)

    nsteps = length(t)          # available steps (including counting
                                # the new n+1'st step)
    h  = t[end]-t[end-1]        # current step size

    if nsteps < k+2
        error("errorEstimates called with too few steps.")
    end

    # estimates the error for orders [k-2,k-1,k]
    errors = zeros(eltype(t),k)

    for i = max(k-2,1):k
        # @todo can this be optimized by inlining the body of
        # interpolateHighestDerivative?
        maxd=interpolateHighestDerivative(t[end-(i+1):end],y[end-(i+1):end])
        errors[i]=h^(i+1)*normy(maxd)
    end

    # compute the estimate the k+1 order only if all the steps were
    # made using the same stepsize
    if nsteps >= k+3
        hn = diff(t[end-(k+2):end])
        if all(hn.==hn[1])
            maxd=interpolateHighestDerivative(t[end-(k+2):end],y[end-(k+2):end])
            push!(errors,h^(k+2)*normy(maxd))
        end
    end

    # return error estimates (this is roughly [ERKM2,ERKM1,ERK,ERKP1]
    # from DASSL)
    return errors

end

# t is an array [t_1,...,t_n] of length n
# y is a matrix [y_1,...,y_n] of size k x l
# h_next is a size of next step
# F encodes the DAE: F(t,y,y')=0
# jd is a bunch of auxilary data saved between steps (jacobian and last coefficient 'a')
# wt is a vector of weights of the norm
function stepper(ord::Integer,
                 t::AbstractVector,
                 y::AbstractVector,
                 dy::AbstractVector,
                 h_next::Real,
                 F,
                 jd::JacData,
                 computejac,
                 wt,
                 normy,
                 maxorder::Integer)

    l        = length(y[1])        # the number of dependent variables

    # @todo this should be the view of the tail of the arrays t and y
    tk = t[end-ord+1:end]
    yk = y[end-ord+1:end]

    t_next   = tk[end]+h_next

    # I think there is an error in the book, the sum should be taken
    # from j=1 to k+1 instead of j=1 to k
    alphas = -sum([1/j for j=1:ord])

    if length(y) == 1
        # this is the first step, we initialize y0 and dy0 with
        # initial data provided by user
        dy0 = dy[1]
        y0  = y[1]+h_next*dy[1]
        # We assume that dy[1] is given, so we can use Hermite
        # interpolation to get a better approximation to the
        # derivative at point y: dy=a*y+b=2(y-y[1])/h-dy[1]
        a=2/h_next
        b=-dy[1]-2*y[1]/h_next
    else
        # we use predictor to obtain the starting point for the
        # modified newton method
        #
        # @todo I should optimize the following functions to return a
        # tuple (y0,dy0)
        dy0 = interpolateDerivativeAt(tk,yk,t_next)
        y0  = interpolateAt(tk,yk,t_next)
        a=-alphas/h_next
        b=dy0-a*y0
    end

    # f_newton is supplied to the modified Newton method.  Zeroes of
    # f_newton give the corrected value of the next step "yc"
    f_newton(yc)=F(t_next,yc,a*yc+b)

    # if called, this function computes the jacobian of f_newton at
    # the point y0 (via first order finite differences or
    # user-supplied functions)
    jac_new()=computejac(t_next,y0,dy0,a)

    # we compute the corrected value "yc", updating the gradient if necessary
    (status,yc,jd)=corrector(jd,       # old coefficient a and jacobian
                             a,        # current coefficient a
                             jac_new,  # this function is called when new jacobian is needed
                             y0,       # starting point for modified newton
                             f_newton, # we want to find zeroes of this function
                             normy)     # the norm used to estimate error needs weights

    alpha = Array{eltype(t)}(undef,ord+1)

    for i = 1:ord
        alpha[i] = h_next/(t_next-t[end-i+1])
    end

    if length(t) >= ord+1
        t0 = t[end-ord]
    elseif length(t) >= 2
        # @todo we choose some arbitrary value of t[0], here t[0]:=t[1]-(t[2]-h[1])
        h1 = t[2]-t[1]
        t0 = t[1]-h1
    else
        t0 = t[1]-h_next
    end

    alpha[ord+1] = h_next/(t_next-t0)

    alpha0 = -sum(alpha[1:ord])
    M      =  max(alpha[ord+1],abs.(alpha[ord+1]+alphas-alpha0))
    err::eltype(t) =  normy((yc-y0))*M

    # status<0 means the modified Newton method did not converge
    # err is the local error estimate from taking the step
    # yc is the estimated value at the next step
    return (status, err, yc, a*yc+b, jd)

end

# returns the corrected value yc and status.  If needed it updates
# the jacobian g_old and a_old.

function corrector(jd::JacData,
                   a_new::T,
                   jac_new,
                   y0::AbstractVector{Ty},
                   f_newton,
                   normy) where {T,Ty}

    # if jd.a == 0 the new jacobian is always computed, independently
    # of the value of a_new
    if abs((jd.a-a_new)/(jd.a+a_new)) > 1/4
        # old jacobian wouldn't give fast enough convergence, we have
        # to compute a current jacobian
        jd = JacData(a_new,jac_new())
        # run the corrector
        (status,yc)=newton_iteration( x->(-(jd.jac\f_newton(x))), y0, normy)
    else
        # old jacobian should give reasonable convergence
        c::T=2*jd.a/(a_new+jd.a) # factor "c" is used to speed up the
                                 # convergence when using an old
                                 # jacobian

        # reusing the old jacobian
        (status,yc)=newton_iteration( x->(-c*(jd.jac\f_newton(x))), y0, normy)

        if status < 0
            # the corrector did not converge, so we recompute jacobian and try again
            jd = JacData(a_new,jac_new())
            # run the corrector again
            (status,yc)=newton_iteration( x->(-(jd.jac\f_newton(x))), y0, normy)
        end
    end

    return (status,yc,jd)

end

# this function iterates f until it finds its fixed point, starting
# from f(y0).  The result either satisfies normy(yn-f(yn))=0+... or is
# set back to y0.  Status tells if the fixed point was obtained
# (status==0) or not (status==-1).
function newton_iteration(f,
                          y0::AbstractVector{T},
                          normy) where T<:Number

    # first guess comes from the predictor method, then we compute the
    # second guess to get the norm1

    delta::typeof(y0)=f(y0)
    norm1::T=normy(delta)
    yn=y0+delta

    # after the first iteration the normy turned out to be very small,
    # terminate and return the first correction step

    ep    = eps(eltype(abs.(y0))) # this is the epsilon for type y0

    if norm1 < 100*ep*normy(y0)
        status=0
        return(status,yn)
    end

    # maximal number of iterations is set by dassl algorithm to 4

    for i=1:MAXIT

        delta=f(yn)
        normn::T=normy(delta)
        rho=(normn/norm1)^(1/i)
        yn=yn+delta

        # iteration failed to converge

        if rho > 9/10
            status=-1
            return(status,y0)
        end

        err=rho/(1-rho)*normn

        # iteration converged successfully

        if err < 1/3
            status=0
            return(status,yn)
        end

    end

    # unable to converge after 4 iterations

    status=-1
    return(status,y0)
end

function dassl_norm(v, wt)
    norm(v./wt)/sqrt(length(v))
end

function dassl_weights(y,reltol,abstol)
    @. reltol*abs(y)+abstol
end

# returns the value of the interpolation polynomial at the point x0
function interpolateAt(x::AbstractVector{T},
                       y::AbstractVector,
                       x0::T) where T<:Real

    if length(x)!=length(y)
        error("x and y have to be of the same size.")
    end

    n = length(x)
    p = zero(y[1])

    for i=1:n
        Li =one(T)
        for j=1:n
            if j==i
                continue
            else
                Li*=(x0-x[j])/(x[i]-x[j])
            end
        end
        p+=Li*y[i]
    end
    return p
end

# returns the value of the derivative of the interpolation polynomial
# at the point x0
function interpolateDerivativeAt(x::AbstractVector{T},
                                 y::AbstractVector,
                                 x0::T) where T<:Real

    if length(x)!=length(y)
        error("x and y have to be of the same size.")
    end

    n = length(x)
    p = zero(y[1])

    for i=1:n
        dLi=zero(T)
        for k=1:n
            if k==i
                continue
            else
                dLi1=one(T)
                for j=1:n
                    if j==k || j==i
                        continue
                    else
                        dLi1*=(x0-x[j])/(x[i]-x[j])
                    end
                end
                dLi+=dLi1/(x[i]-x[k])
            end
        end
        p+=dLi*y[i]
    end
    return p
end

# if the interpolating polynomial is given as
# p(x)=a_{k-1}*x^{k-1}+...a_1*x+a_0 then this function returns the
# k-th derivative of p, i.e. (k-1)!*a_{k-1}
function interpolateHighestDerivative(x::AbstractVector,
                                      y::AbstractVector)

    if length(x)!=length(y)
        error("x and y have to be of the same size.")
    end

    n = length(x)
    p = zero(y[1])

    for i=1:n
        Li =one(eltype(x))
        for j=1:n
            if j==i
                continue
            else
                Li*=1/(x[i]-x[j])
            end
        end
        p+=Li*y[i]
    end
    return factorial(n-1)*p
end

# generate a function that computes approximate jacobian using forward
# finite differences
function numerical_jacobian(F,reltol,abstol,weights)
    function numjac(t, y, dy, a)
        ep      = eps(one(t))   # this is the machine epsilon
        h       = 1/a           # h ~ 1/a
        wt      = weights(y,reltol,abstol)
        # delta for approximation of jacobian.  I removed the
        # sign(h_next*dy0) from the definition of delta because it was
        # causing trouble when dy0==0 (which happens for ord==1)
        edelta  = diagm(0=>max.(abs.(y),abs.(h*dy),wt)*sqrt(ep))

        b=dy-a*y
        f(y1) = F(t,y1,a*y1+b)

        n   = length(y)
        jac = Array{eltype(y)}(undef,n,n)
        for i=1:n
            jac[:,i]=(f(y+edelta[:,i])-f(y))/edelta[i,i]
        end

        jac
    end
end

end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(halfquadruped)
v1 = zeros(halfquadruped.nq)

# ## Time
h = 0.01
T = 100

# ## Simulator
s = Simulator(halfquadruped, T, h=h)

# ## Simulate
simulate!(s, q1, v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)


text:
using QContinuousMeasurement
using Test
using Random
using QContinuousMeasurement
using Test
using Random

@testset "collective" begin
    @testset "Nj = 1" begin
        # We check that Local and Collective dephasing are the same
        # for a single spin.

        params = ModelParameters(Nj=1, kind=1, kcoll=0, Gamma=1)
        model = LocalDephasingModel(params)
        state = blockdiag_css(params.Nj)

        params_coll = ModelParameters(Nj=1, kind=0, kcoll=2, Gamma=1)
        model_coll = CollectiveDephasingModel(params_coll)
        state_coll = fixedj_css(params_coll.Nj)

        seed = 1

        Random.seed!(seed)
        res = simulate_trajectory(model, state)

        Random.seed!(seed)
        res_coll = simulate_trajectory(model_coll, state_coll)

        for k in keys(res)
            @testset "$k" begin
                @test getfield(res, k) ≈ getfield(res_coll, k)
            end
        end
    end

    @testset "Nj = 5" begin
        # We check that <Jx> is the same for local and collective
        # when there is no monitoring

        params = ModelParameters(Nj=5, kind=1, kcoll=0, Gamma=0)
        model = LocalDephasingModel(params)
        state = blockdiag_css(params.Nj)

        params_coll = ModelParameters(Nj=5, kind=0, kcoll=2, Gamma=0)
        model_coll = CollectiveDephasingModel(params_coll)
        state_coll = fixedj_css(params_coll.Nj)

        seed = 1

        Random.seed!(seed)
        res = simulate_trajectory(model, state)

        Random.seed!(seed)
        res_coll = simulate_trajectory(model_coll, state_coll)

        @testset "Jx" begin
            @test getfield(res, :Jx) ≈ getfield(res_coll, :Jx) rtol = 10 * params.dt
        end
    end
end


text:
# The prime factors of 13195 are 5, 7, 13 and 29.
#
# What is the largest prime factor of the number 600851475143?

using ProjectEulerSolutions

# Call the prime factors function, return the last.
function p003solution(n::Integer=3)::Integer
    return primefactors(n)[end]
end

p003 = Problems.Problem(p003solution)

Problems.benchmark(p003, 600_851_475_143)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project