text:
# If a library is not found use Pkg.add("Datetime") from the Julia command line to add
# to the libraries available on this machine.
# The very first time you will get Unable to read directory METADATA
# so set the directory ENV["JULIA_PKGDIR"] = "c:/Users/simag001/.julia"
# then Pkg.Init()
# Use the Datetime library for date and time manipulations
# This is deprecated, should be using Dates and Julia 0.3
# Check version using versioninfo()
using Dates
using DataArrays, DataFrames
using GLM

include("temp_func.jl")

function predict_with_offset(temp_offset::Float64,min_month::Int,max_month::Int,airport_code::ASCIIString)
# Run through the whole five years and calculate the (1-prob0) value, ie expected value
currdate=DateTime(2009,5,1,0,0,0,0)
expected=0.0
while currdate < DateTime(2014,5,1,0,0,0,0)
  if (month(currdate)>=min_month && month(currdate)<=max_month)
  # for reach currdate calculate the expected value
  beta=coef(MNL)[1]
  beta+=coef(MNL)[2]*(year(currdate)-2005)
  # Handle the hour curve
  hourfrac=hour(currdate)+minute(currdate)/60
  beta+=coef(MNL)[3]*(min(1-0,max(0,hourfrac-0)))
  beta+=coef(MNL)[4]*(min(2-1,max(0,hourfrac-1)))
  beta+=coef(MNL)[5]*(min(4-3,max(0,hourfrac-3)))
  beta+=coef(MNL)[6]*(min(5-4,max(0,hourfrac-4)))
  beta+=coef(MNL)[7]*(min(6-5,max(0,hourfrac-5)))
  beta+=coef(MNL)[8]*(min(7-6,max(0,hourfrac-6)))
  beta+=coef(MNL)[9]*(min(8-7,max(0,hourfrac-7)))
  beta+=coef(MNL)[10]*(min(9-8,max(0,hourfrac-8)))
  beta+=coef(MNL)[11]*(min(10-9,max(0,hourfrac-9)))
  beta+=coef(MNL)[12]*(min(11-10,max(0,hourfrac-10)))
  beta+=coef(MNL)[13]*(min(12-11,max(0,hourfrac-11)))
  beta+=coef(MNL)[14]*(min(13-12,max(0,hourfrac-12)))
  beta+=coef(MNL)[15]*(min(16-15,max(0,hourfrac-15)))
  beta+=coef(MNL)[16]*(min(17-16,max(0,hourfrac-16)))
  beta+=coef(MNL)[17]*(min(18-17,max(0,hourfrac-17)))
  beta+=coef(MNL)[18]*(min(19-18,max(0,hourfrac-18)))
  beta+=coef(MNL)[19]*(min(20-19,max(0,hourfrac-19)))
  beta+=coef(MNL)[20]*(min(21-20,max(0,hourfrac-20)))
  beta+=coef(MNL)[21]*(min(22-21,max(0,hourfrac-21)))
  beta+=coef(MNL)[22]*(min(23-22,max(0,hourfrac-22)))
  beta+=coef(MNL)[23]*(min(24-23,max(0,hourfrac-23)))
  # section for Date variables
  if dayofweek(currdate)==Dates.Monday
    beta+=coef(MNL)[24]
  end
  if dayofweek(currdate)==Dates.Tuesday
    beta+=coef(MNL)[25]
  end
  if dayofweek(currdate)==Dates.Wednesday
    beta+=coef(MNL)[26]
  end
  if dayofweek(currdate)==Dates.Thursday
    beta+=coef(MNL)[27]
  end
  if dayofweek(currdate)==Dates.Friday
    beta+=coef(MNL)[28]
  end
  if dayofweek(currdate)==Dates.Saturday
    beta+=coef(MNL)[29]
  end
  # temperature section
  if airport_code=="ADL"
    currtemp=temp_offset+get_adelaide_temp(currdate)
  else
    currtemp=temp_offset+get_tokyo_temp(currdate)
  end
  beta+=coef(MNL)[30]*min(15-12,max(0,currtemp-12))
  beta+=coef(MNL)[31]*min(18-15,max(0,currtemp-15))
  beta+=coef(MNL)[32]*min(22-18,max(0,currtemp-18))
  beta+=coef(MNL)[33]*min(28-25,max(0,currtemp-25))
  beta+=coef(MNL)[34]*min(31-28,max(0,currtemp-28))
  beta+=coef(MNL)[35]*min(37-34,max(0,currtemp-34))
  beta+=coef(MNL)[36]*min(40-37,max(0,currtemp-37))
  # dew point
  if airport_code=="ADL"
    beta+=coef(MNL)[37]*get_adelaide_dewp(currdate)
  else
    beta+=coef(MNL)[37]*get_tokyo_dewp(currdate)
  end
  # wind speed
  if airport_code=="ADL"
    beta+=coef(MNL)[38]*get_adelaide_wind_speed(currdate)
  else
    beta+=coef(MNL)[38]*get_tokyo_wind_speed(currdate)
  end
  # for Poisson distribution we expected value is equal to lambda
  expected+=exp(beta)
  end
  currdate=currdate+Dates.Minute(1)
  end
  return expected
end

source_file="mfs.csv"
# Make sure we can read our data file
if (!isfile(source_file))
   # We need to fetch the file
   download("http://data.sa.gov.au/storage/f/2014-06-26T05%3A25%3A56.821Z/samfs-data-may-01-2009-to-30-apr-2014.csv", source_file)
   if (!isfile(source_file))
      error("Tried to download $source_file and that didn't work, giving up")
   end
end

# Read in the file
d2=readlines(open(source_file));
# Create a dictionary to hold all incident times
calllog=Dict{DateTime,Int}()
# Iterate through the file loading each elem
incidents_loaded=0
for thiselem in d2[2:end]
 eachelem=split(chomp(thiselem),',',4)
 if eachelem[4]!="MUTUAL AID GIVEN/CHANGE OF QUARTERS"
   thisdate=DateTime(int(eachelem[2][1:4]),int(eachelem[2][6:7]),int(eachelem[2][9:10]),int(eachelem[2][12:13]),int(eachelem[2][15:16]),int(eachelem[2][18:19]),0)
    if ! haskey(calllog, thisdate)
      calllog[thisdate]=1
    else
      calllog[thisdate]+=1
    end
   incidents_loaded+=1
  end
 end
 # Iterate through each time loading a fire or not
thisdate=DateTime(2009,5,1,0,0,0,0)
zeroskipped=0
linesprinted=0
FireVec=Int[]
DateVec=DateTime[]
TempVec=Float64[]
DewpVec=Float64[]
WindSpeedVec=Float64[]
while thisdate < DateTime(2014,5,1,0,0,0,0)
  zeroskipped=zeroskipped+1
  # 12 15 18 22
  current_temp=get_adelaide_temp(thisdate)
  current_dewp=get_adelaide_dewp(thisdate)
  current_wind_speed=get_adelaide_wind_speed(thisdate)
  if (haskey(calllog,thisdate))
    if linesprinted < 1
      # use the number of incidents in this time slot
      push!(FireVec,calllog[thisdate])
      push!(DateVec,thisdate)
      push!(TempVec,current_temp)
      push!(DewpVec,current_dewp)
      push!(WindSpeedVec, current_wind_speed)
      linesprinted=linesprinted+1
    end
   end
  if zeroskipped >=1
    if linesprinted < 1
        push!(FireVec,0)
        push!(DateVec,thisdate)
        push!(TempVec,current_temp)
        push!(DewpVec,current_dewp)
        push!(WindSpeedVec, current_wind_speed)
      end
      zeroskipped=0
      linesprinted=0
    end
 thisdate=thisdate+Dates.Minute(1)
end
df= DataFrame()
df[:Fire]=FireVec
df[:Incdate]=DateVec
df[:dayofweek]=dayofweek(DateVec)
df[:year]=year(DateVec)-2005
df[:hour]=hour(DateVec)
hourfrac=hour(DateVec)+minute(DateVec)/60
df[:hourfrac]=hourfrac
df[:hourfrac2]=hourfrac.^2
df[:hourfrac3]=hourfrac.^3
df[:hourfrac4]=hourfrac.^4
df[:hourfrac5]=hourfrac.^5
df[:hourfrac6]=hourfrac.^6
df[:hourfrac7]=hourfrac.^7
df[:hourfrac8]=hourfrac.^8
df[:hourfrac9]=hourfrac.^9
df[:hourfrac10]=hourfrac.^10
df[:hourfrac11]=hourfrac.^11
df[:hour0to1]=min(1-0,max(0,hourfrac-0))
df[:hour1to2]=min(2-1,max(0,hourfrac-1))
df[:hour2to3]=min(3-2,max(0,hourfrac-2))
df[:hour3to4]=min(4-3,max(0,hourfrac-3))
df[:hour4to5]=min(5-4,max(0,hourfrac-4))
df[:hour5to6]=min(6-5,max(0,hourfrac-5))
df[:hour6to7]=min(7-6,max(0,hourfrac-6))
df[:hour7to8]=min(8-7,max(0,hourfrac-7))
df[:hour8to9]=min(9-8,max(0,hourfrac-8))
df[:hour9to10]=min(10-9,max(0,hourfrac-9))
df[:hour10to11]=min(11-10,max(0,hourfrac-10))
df[:hour11to12]=min(12-11,max(0,hourfrac-11))
df[:hour11to13]=min(13-11,max(0,hourfrac-11))
df[:hour12to13]=min(13-12,max(0,hourfrac-12))
df[:hour13to14]=min(14-13,max(0,hourfrac-13))
df[:hour14to15]=min(15-14,max(0,hourfrac-14))
df[:hour15to16]=min(16-15,max(0,hourfrac-15))
df[:hour16to17]=min(17-16,max(0,hourfrac-16))
df[:hour17to18]=min(18-17,max(0,hourfrac-17))
df[:hour18to19]=min(19-18,max(0,hourfrac-18))
df[:hour19to20]=min(20-19,max(0,hourfrac-19))
df[:hour20to21]=min(21-20,max(0,hourfrac-20))
df[:hour21to22]=min(22-21,max(0,hourfrac-21))
df[:hour22to23]=min(23-22,max(0,hourfrac-22))
df[:hour23to24]=min(24-23,max(0,hourfrac-23))
df[:month]=month(DateVec)
df[:temperature]=TempVec
df[:temp2]=TempVec.^2
df[:temp3]=TempVec.^3
df[:temp4]=TempVec.^4
df[:temp5]=TempVec.^5
df[:temp6]=TempVec.^6
df[:temp7]=TempVec.^7
df[:temp8]=TempVec.^8
df[:temp9]=TempVec.^9
df[:temp10]=TempVec.^10
df[:IsMonday]=0+(df[:dayofweek].==1)
df[:IsTuesday]=0+(df[:dayofweek].==2)
df[:IsWednesday]=0+(df[:dayofweek].==3)
df[:IsThursday]=0+(df[:dayofweek].==4)
df[:IsFriday]=0+(df[:dayofweek].==5)
df[:IsSaturday]=0+(df[:dayofweek].==6)
df[:IsSunday]=0+(df[:dayofweek].==7)
df[:tempband1]=max(0,12-TempVec)
df[:tempband2]=min(15-12,max(0,TempVec-12))
df[:tempband3]=min(18-15,max(0,TempVec-15))
df[:tempband4]=min(22-18,max(0,TempVec-18))
df[:tempband5a]=min(25-22,max(0,TempVec-22))
df[:tempband5b1]=min(28-25,max(0,TempVec-25))
df[:tempband5b2]=min(31-28,max(0,TempVec-28))
df[:tempband5b3]=min(34-31,max(0,TempVec-31))
df[:tempband5b4]=min(37-34,max(0,TempVec-34))
df[:tempband5b5]=min(40-37,max(0,TempVec-37))
df[:tempband5bz]=min(100-40,max(0,TempVec-40))
df[:dew_point]=DewpVec
df[:dew_point2]=DewpVec.^2
df[:relative_humidity]=100-5(TempVec-DewpVec)
df[:wind_speed]=WindSpeedVec
df[:wind_speed2]=WindSpeedVec.^2
# writetable("comb.csv",df)
# We know we are finished and how many incidents there were
println("$incidents_loaded incidents loaded")
MNL= glm(Fire ~ year + hourfrac+hourfrac2+hourfrac3+hourfrac4+hourfrac5+hourfrac6+hourfrac7+hourfrac8+hourfrac9 +hourfrac10+ IsMonday + IsTuesday + IsWednesday+ IsThursday + IsFriday + IsSaturday +temperature+temp2+temp3+temp4+temp5+temp6 +dew_point+dew_point2+wind_speed+wind_speed2 , df, Poisson())
# Now go through and predict manually to extract the expected number of fires with a change in temperature
#println("Expected Value for 0 ", predict_with_offset(0.0,1,12,"ADL"))
#println("Expected Value for 1 ", predict_with_offset(1.0,1,12,"ADL"))
#println("Expected Value for 2 ", predict_with_offset(2.0,1,12,"ADL"))
#println("Expected Value for -1 ", predict_with_offset(-1.0,1,12,"ADL"))
#println("Expected Value for 0 ", predict_with_offset(0.0,1,12,"HND"))
#println("Expected Value for 1 ", predict_with_offset(1.0,1,12,"HND"))
#println("Expected Value for 2 ", predict_with_offset(2.0,1,12,"HND"))
#println("Expected Value for -1 ", predict_with_offset(-1.0,1,12,"HND"))

max_stars_repo_name: softgrow/mfsmodel-weather

max_issues_repo_name: softgrow/mfsmodel-weather

max_forks_repo_name: softgrow/mfsmodel-weather


text:
function epochsec2datetime(epochsec::Int64)::DateTime
    return Dates.epochms2datetime(Int64(epochsec * 1000))
end

function plot_schedule(
    scheduled_jobs::Vector{Job},
    machines::Vector{Machine};
    xlims = (-Inf, Inf),
    highlight_job_id = nothing,
    size = (600, 400),
)

    # axes
    plot(
        repeat([Inf], inner = length(machines)),
        reverse([m.name for m in machines]),
        yticks = (0.5:1:(length(machines)-0.5), reverse([m.name for m in machines])),
        ylims = (0, length(machines)),
        size = size,
    )

    for job in scheduled_jobs
        # operations
        color = job.id
        if !isnothing(highlight_job_id)
            if job.id == highlight_job_id
                color = :darkblue
            else
                color = :lightgray
            end
        end
        for op in job.operations
            xs = [epochsec2datetime(op.S), epochsec2datetime(end_time(op))]
            ys = repeat([machines[op.E].name], inner = 2)
            plot!(xs, ys, linewidth = 4, linecolor = color)
        end

        # operation dependencies
        for edge in job.dag
            xs = [epochsec2datetime(end_time(job.operations[edge[1]])), epochsec2datetime(job.operations[edge[2]].S)]
            ys = [machines[job.operations[edge[1]].E].name, machines[job.operations[edge[2]].E].name]
            plot!(xs, ys, linestyle = :dash, linecolor = :gray)
        end
    end

    plot!(
        legend = :none,
        title = "Scheduling Result",
        xlabel = "Date/Time",
        ylabel = "Machine",
        xlims = xlims,
    )
end

ext: jl


text:
import .Utils: hill, hillr, exprel, mM, μM, ms, minute, Hz, mV, iVT, iVmtx, iCmt, iVi, iVimtx, F_M

@variables t

#=
Glucokinase (GK) Reactions:

- Glc + ATP => G6P (+ ADP)
- G6P + ATP => FBP (+ ADP) =>> 2G3P (+ ADP)
=#
@variables Glc(t) ATP_c(t) J_GK(t)
@parameters VmaxGK = 0.011mM * Hz KatpGK = 0.5mM KglcGK = 7mM NGK = 1.7 GlcConst = 5mM

#=
Glyceraldehydes 3-phosphate dehydrogenase (GPD), lumped with pyruvate kinase (PK)

G3P + NAD + 2ADP => Pyr + NADH + 2ATP
=#
@variables J_GPD(t) G3P(t) NAD_c(t) NADH_c(t) ADP_c(t)
@parameters VmaxGPD = 0.5mM * Hz Kg3pGPD = 0.2mM KnadGPD = 0.09 KadpGPD = 0.020mM

#=
Lactate dehydrogenase (LDH)
The rate of lactate output is approximately 5% of the glycolytic flux at 8 mM glucose.

Pyr + NADH => Lactate + NAD
=#
@variables J_LDH(t) Pyr(t)
@parameters VmaxLDH = 1.2mM * Hz KpyrLDH = 47.5μM KnadhLDH = 1.0

# Cytosolic calcium
@variables Ca_c(t)
@parameters RestingCa = 0.09μM ActivatedCa = 0.25μM NCac = 4 KatpCac = 25.0

#=
Adenylate kinase:

2ADP <=> ATP + AMP
=#
@variables J_ADK(t) AMP_c(t)
@parameters kfAK = 1000Hz / mM kEqAK = 0.931

# Activity of AMPK
@variables AMPKactivity(t)
@parameters kAMPK = 0.07

#=
Pyruvaye dehydrogenase (PDH) as well as electron transport chain (ETC)

Pyr + 4.6NAD => CO2 + 4.6NADH
=#
@variables J_PDH(t) J_DH(t) J_ETC(t) NAD_m(t) NADH_m(t) Ca_m(t) rPDH(t)
@parameters VmaxPDH = 300μM * Hz KpyrPDH = 47.5μM KnadPDH = 81.0 U1PDH = 1.5 U2PDH = 1.1 KcaPDH = 0.05μM

function j_pdh(pyr, nad_m, nadh_m, ca_m, VMAX, K_PYR, K_NAD, U1, U2, K_CA)
    c = (hillr(ca_m, K_CA))^2
    fpCa = hillr(U2 * (1 + U1 * c))
    fNAD = hill(nad_m, nadh_m * K_NAD)
    fPyr = hill(pyr, K_PYR)
    jPDH = VMAX * rPDH * fPyr * fNAD * fpCa
    return jPDH
end

#=
Electron trasnport chain (ETC)

NADH + O2 + 10H(mito) => NAD + 10H(cyto)
=#
@variables J_HR(t) J_O2(t) ΔΨm(t) rETC(t)
@parameters VmaxETC = 22mM * Hz KnadhETC = 3mM KaETC = -4.92E-3 / mV KbETC = -4.43E-3 / mV

# Proton leak
@variables J_HL(t) rHL(t)
@parameters pHleak = 2.4μM * Hz kvHleak = 0.0305 / mV

#=
F1Fo ATPase (ATP synthase) plus ANT

3ADP + 8H(cyto) => 3ATP + 8H(mito)
=#
@variables J_HF(t) J_ANT(t) rF1(t)
@parameters VmaxF1 = 8mM * Hz KadpF1 = 20μM KvF1 = 131.4mV KcaF1 = 0.165μM FmgadpF1 = 0.055

#=
Mitochondrial calcium uniporter (MCU)

Ca(cyto) <=> Ca(mito)
=#
@variables J_MCU(t)
@parameters PcaMCU = 4Hz
function j_uni(ca_m, ca_c, ΔΨ, P_CA)
    zvfrt = 2 * iVT * ΔΨ
    em1 = expm1(zvfrt)
    return P_CA * exprel(zvfrt, em1) * (0.341 * ca_c * (1 + em1) - 0.200 * ca_m)
end

#=
Mitochondrial sodium calcium exchanger (NCLX)

Ca(mito) + 2Na(cyto) <=> Ca(cyto) + 2Na(mito)
=#
@variables J_NCLX(t)
@parameters Na_c = 10mM Na_m = 5mM VmaxNCLX = 75μM * Hz KnaNCLX = 8.2mM KcaNCLX = 8μM
function j_nclx(ca_m, ca_c, na_m, na_c, VMAX, K_CA, K_NA)
    A = (na_c / K_NA)^2
    P = (na_m / K_NA)^2
    B = ca_m / K_CA
    Q = ca_c / K_CA
    AB = A * B
    PQ = P * Q
    return VMAX * (AB - PQ) / (1 + A + B + P + Q + AB + PQ)
end

#=
NADH shuttle

NADH(cyto) + NAD(mito) <=> NADH(mito) + NAD(cyto)
=#
@variables J_NADHT(t)
@parameters VmaxNADHT = 50μM * Hz Ktn_c = 0.002 Ktn_m = 16.78

# Fission-fusion rates
@variables x[1:3](t) x13r(t) degavg(t) v[1:2](t)
@parameters Kfiss1 = inv(10minute) Kfuse1 = Kfiss1 Kfiss2 = 1.5Kfiss1 Kfuse2 = 0.5Kfuse1

# Baseline consumption rates
@parameters kNADHc = 0.1Hz kNADHm = 0.1Hz kATP = 0.04Hz kATPCa = 90Hz / mM kG3P = 0.01Hz kPyr = 0.01Hz

# Conservation relationships
@parameters ΣAc = 4.5mM Σn_c = 2mM Σn_m = 2.2mM

const DEFAULT_U0 = Dict(
    G3P => 2.8μM,
    Pyr => 8.5μM,
    NADH_c => 1μM,
    NADH_m => 60μM,
    ATP_c => 4000μM,
    AMP_c => 0μM,
    Ca_m => 0.250μM,
    ΔΨm => 100mV,
    x[2] => 0.25,
    x[3] => 0.05
)

function make_model(;
    name,
    simplify=true,
    cacrhs=RestingCa + ActivatedCa * hill(ATP_c / KatpCac, ADP_c, NCac),
    glcrhs=GlcConst,
    rpdh=1,
    retc=1,
    rf1=1,
    rhleak=1
)
    D = Differential(t)
    eqs = [
        # Reactions
        J_GK ~ VmaxGK * hill(ATP_c, KatpGK) * hill(Glc, KglcGK),
        J_GPD ~ VmaxGPD * hill(ADP_c, KadpGPD) * hill(NAD_c / KnadGPD, NADH_c) * hill(G3P, Kg3pGPD),
        J_LDH ~ VmaxLDH * hill(Pyr, KpyrLDH) * hill(NADH_c / KnadhLDH, NAD_c),
        J_ADK ~ kfAK * (ADP_c * ADP_c - ATP_c * AMP_c / kEqAK),
        AMPKactivity ~ hill(AMP_c / ATP_c, kAMPK),
        J_PDH ~ j_pdh(Pyr, NAD_m, NADH_m, Ca_m, VmaxPDH, KpyrPDH, KnadPDH, U1PDH, U2PDH, KcaPDH),
        J_ETC ~ J_PDH,
        J_DH ~ 4.6 * J_ETC,
        J_HR ~ VmaxETC * rETC * hill(NADH_m, KnadhETC) * (1 + KaETC * ΔΨm) / (1 + KbETC * ΔΨm),
        J_O2 ~ 0.1 * J_HR,
        J_HL ~ pHleak * rhleak * exp(kvHleak * ΔΨm),
        J_HF ~ VmaxF1 * rF1 * hill(FmgadpF1 * ADP_c, KadpF1, 2) * hill(ΔΨm, KvF1, 8) * (1 - exp(-Ca_m / KcaF1)),
        J_ANT ~ J_HF / 3,
        J_MCU ~ j_uni(Ca_m, Ca_c, ΔΨm, PcaMCU),
        J_NCLX ~ j_nclx(Ca_m, Ca_c, Na_m, Na_c, VmaxNCLX, KcaNCLX, KnaNCLX),
        J_NADHT ~ VmaxNADHT * hill(NADH_c / Ktn_c, NAD_c) * hill(NAD_m / Ktn_m, NADH_m),
        v[1] ~ Kfuse1 * J_ANT / J_HL * x[1] * x[1] - Kfiss1 * x[2],
        v[2] ~ Kfuse2 * J_ANT / J_HL * x[1] * x[2] - Kfiss2 * x[3],
        Glc ~ glcrhs,
        Ca_c ~ cacrhs,
        rPDH ~ rpdh,
        rETC ~ retc,
        rF1 ~ rf1,
        rHL ~ rhleak,
        # Conservation relationships
        ΣAc ~ ATP_c + ADP_c + AMP_c,
        Σn_c ~ NADH_c + NAD_c,
        Σn_m ~ NADH_m + NAD_m,
        1 ~ x[1] + 2x[2] + 3x[3],
        # Observables
        x13r ~ x[1] / x[3],
        degavg ~ (x[1] + 2x[2] + 3x[3]) / (x[1] + x[2] + x[3]),
        # State variables
        D(NADH_m) ~ iVmtx * (J_DH + J_NADHT - J_O2) - kNADHm * NADH_m,
        # D(NAD_m) ~ # Conserved
        D(NADH_c) ~ iVi * (J_GPD - J_NADHT - J_LDH) - kNADHc * NADH_c,
        # D(NAD_c) ~ # Conserved
        D(Ca_m) ~ iVmtx * F_M * (J_MCU - J_NCLX),
        D(ΔΨm) ~ iCmt * (J_HR - J_HF - J_HL - J_ANT - 2J_MCU),
        D(G3P) ~ iVi * (2J_GK - J_GPD) - kG3P * G3P,
        D(Pyr) ~ iVimtx * (J_GPD - J_PDH - J_LDH) - kPyr * Pyr,
        D(ATP_c) ~ iVi * (-2J_GK + 2J_GPD + J_ANT + J_ADK) - ATP_c * (kATP + kATPCa * Ca_c),
        D(AMP_c) ~ iVi * J_ADK,
        # D(ADP_c) ~ # Conserved
        # D(x[1]) ~ # Conserved
        D(x[2]) ~ v[1] - v[2],
        D(x[3]) ~ v[2],
    ]
    sys = ODESystem(eqs, t; name, defaults=DEFAULT_U0)

    if simplify
        sys = structural_simplify(sys)
    end

    return sys
end

max_stars_repo_name: ntumitolab/mitodyn-ode

max_issues_repo_name: ntumitolab/mitodyn-ode

max_forks_repo_name: NTUMitoLab/mitodyn-ode


text:
mutable struct KKTErrorEvaluator{G,C,J}
    f::Vector{Float64}
    c::Vector{Float64}

    jac::SparseMatrixCSC{Float64,Int}
    jac_view::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}

    grad!::G
    con!::C
    jac!::J
end

function KKTErrorEvaluator(m::Model)

    obj = m.obj
    con = m.con
    p = m.p

    grad = Gradient(obj)
    jacobian,jac_sparsity = SparseJacobian(con)

    _obj = @inline x->non_caching_eval(obj,x,p)
    _con! = @inline (y,x)->non_caching_eval(con,y,x,p)
    _grad! = @inline function (y,x)
        y .= 0
        obj(x,p)
        non_caching_eval(grad,y,x,p)
    end
    _jac! = @inline function (J,x)
        J .= 0
        con(DUMMY,x,p)
        non_caching_eval(jacobian,J,x,p)
    end
    _jac_fill_sparsity! = @inline (I,J)->fill_sparsity!(I,J,jac_sparsity)

    nnz_jac = length(jac_sparsity)

    I = Vector{Int}(undef,nnz_jac)
    J = Vector{Int}(undef,nnz_jac)
    _jac_fill_sparsity!(I,J)

    f = similar(m.x)
    c = similar(m.l)

    tmp = sparse(J,I,collect(1:nnz_jac),num_variables(m),num_constraints(m))
    jac = SparseMatrixCSC{Float64,Int64}(tmp.m,tmp.n,tmp.colptr,tmp.rowval,Vector{Float64}(undef,nnz_jac))

    idx = Vector{Int}(undef,length(tmp.nzval))
    for i=eachindex(tmp.nzval)
        idx[tmp.nzval[i]] = i
    end
    jac_view = view(jac.nzval,idx)

    return KKTErrorEvaluator(f,c,jac,jac_view,_grad!,_con!,_jac!)
end

function (evaluator::KKTErrorEvaluator{G,C,J})(x,l,gl) where {G,C,J}

    evaluator.grad!(evaluator.f,x)
    evaluator.con!(evaluator.c,x)
    evaluator.jac!(evaluator.jac_view,x)

    mul!(evaluator.f,evaluator.jac,l,1.,1.)
    evaluator.c .-= gl

    return max(norm(evaluator.f,Inf),norm(evaluator.c,Inf))
end


text:
module DukascopyTicksReader

    export get, DukascopyTicks, CacheDirectory
    import Base: download

    include("utils.jl")
    include("reader.jl")

    using Formatting

    SKIP_ERROR_DEFAULT = false

    const Ticker = AbstractString
    const Tickers = Vector{String}

    abstract type AbstractCache end

    struct NoCache <: AbstractCache
    end

    struct CacheDirectory <: AbstractCache
        dir::AbstractString

        CacheDirectory() = new("")
        CacheDirectory(dir) = new(dir)
    end

    abstract type AbstractDataCacheSource end

    struct DataFromNetwork <: AbstractDataCacheSource
    end

    struct DataFromCache <: AbstractDataCacheSource
    end

    abstract type AbstractDataReader end

    struct DukascopyTicks <: AbstractDataReader
        cache::AbstractCache

        function DukascopyTicks()
            new(CacheDirectory(""))
        end
    end

    function _cache_dir(dr::DukascopyTicks, cache::CacheDirectory, ticker::Ticker, dt::DateTime)
        if cache.dir == ""
            d = Date(dt)
            #dt_round = DateTime(Dates.year(dt), Dates.month(dt), Dates.day(dt), Dates.hour(dt))
            #joinpath(homedir(), "data", "dukascopy", "ticks",
            #               string(Dates.year(d)), string(d)[1:end-3], string(d), ticker)
            joinpath(homedir(), "data", "dukascopy", "ticks",
                       string(Dates.year(d)), string(d)[1:end-3], string(d), Dates.format(dt, "yyyy-mm-dd_HH0000"))
        else
            cache.dir
        end
    end

    function _url(dr::DukascopyTicks, ticker::Ticker, dt::DateTime)
        yy = Dates.year(dt)
        mm = Dates.month(dt)
        dd = Dates.day(dt)
        hh = Dates.hour(dt)
        format("http://www.dukascopy.com/datafeed/{1}/{2:04d}/{3:02d}/{4:02d}/{5:02d}h_ticks.bi5", ticker, yy, mm, dd, hh)
    end

    function _destination_filename(dr::DukascopyTicks, ticker::Ticker, dt::DateTime)
        #hh = Dates.hour(dt)
        #format("{1:02d}h_ticks.bi5", hh)
        ticker * ".bi5"
    end

    function _cache_filename(dr::DukascopyTicks, cache::CacheDirectory, ticker::Ticker, dt::DateTime)
        joinpath(_cache_dir(dr, cache, ticker, dt), _destination_filename(dr, ticker, dt))
    end

    function get(dr::DukascopyTicks, ticker::Ticker, dt::DateTime, ::DataFromCache, cache::CacheDirectory, skip_error::Bool)
        filename = _cache_filename(dr, cache, ticker, dt)
        println("get $ticker for $dt from fname=$filename")
    end

    function get(dr::DukascopyTicks, ticker::Ticker, dt::DateTime, ::DataFromNetwork, cache::CacheDirectory, skip_error::Bool)
        url = _url(dr, ticker, dt)
        filename = _cache_filename(dr, cache, ticker, dt)
        if !ispath(filename)
            mkpath(_cache_dir(dr, cache, ticker, dt))
        end
        println("download $url")
        try
            download(url, filename)
        catch
            if skip_error
                println("skip_error")
            else
                rethrow()
            end
        end
        println("save $ticker for $dt to cache $filename")
    end

    function _is_in_cache(dr::DukascopyTicks, ticker::Ticker, dt::DateTime, cache::CacheDirectory)
        filename = joinpath(_cache_dir(dr, cache, ticker, dt), _destination_filename(dr, ticker, dt))
        isfile(filename)
    end

    function download(dr::DukascopyTicks, ticker::Ticker, dt::DateTime; skip_error=SKIP_ERROR_DEFAULT)
        cache = dr.cache
        if _is_in_cache(dr, ticker, dt, cache)
            get(dr, ticker, dt, DataFromCache(), cache, skip_error)
        else
            get(dr, ticker, dt, DataFromNetwork(), cache, skip_error)
        end
    end

    function get(dr::DukascopyTicks, ticker::Ticker, dt::DateTime; skip_error=SKIP_ERROR_DEFAULT)
        cache = dr.cache
        download(dr, ticker, dt, skip_error=skip_error)
        filename = _cache_filename(dr, cache, ticker, dt)
        reader = TickReader(dt::DateTime, ticker, filename)
        reader
    end

    function download(dr::AbstractDataReader, ticker::Union{Ticker, Tickers}, dt_range::StepRange; skip_error=SKIP_ERROR_DEFAULT)
        for dt in dt_range
            download(dr, ticker, dt, skip_error=skip_error)
        end
    end

    function download(dr::DukascopyTicks, ticker::Union{Ticker, Tickers}, start::DateTime, stop::DateTime; skip_error=SKIP_ERROR_DEFAULT)
        step = Dates.Hour(1)
        download(dr, ticker, start:step:stop-step, skip_error=skip_error)
    end

    function download(dr::AbstractDataReader, tickers::Tickers, dt::DateTime; skip_error=SKIP_ERROR_DEFAULT)
        for ticker in tickers
            download(dr, ticker, dt, skip_error=skip_error)
        end
    end

end # module


text:
"""
Simulates the neoclassical growth model for a given set of solution
coefficients. It simulates for `capT` periods and discards first
`nburn` observations.
"""
function simulate(ncgm::NeoclassicalGrowth, vp::ValueCoeffs,
                  shocks::Vector{Float64}; capT::Int=10_000,
                  nburn::Int=200)
    # Unpack parameters
    kp = 0.0  # Policy holder
    temp = Array{Float64}(undef, n_complete(2, vp.d))

    # Allocate space for k and z
    ksim = Array{Float64}(undef, capT+nburn)
    zsim = Array{Float64}(undef, capT+nburn)

    # Initialize both k and z at 1
    ksim[1] = 1.0
    zsim[1] = 1.0

    # Simulate
    temp = Array{Float64}(undef, n_complete(2, vp.d))
    for t in 2:capT+nburn
        # Evaluate k_t given yesterday's (k_{t-1}, z_{t-1})
        kp = env_condition_kp!(temp, ncgm, vp, ksim[t-1], zsim[t-1])

        # Draw new z and update k using policy above
        zsim[t] = zsim[t-1]^ncgm.ρ * exp(ncgm.σ*shocks[t])
        ksim[t] = kp
    end

    return ksim[nburn+1:end], zsim[nburn+1:end]
end

function simulate(ncgm::NeoclassicalGrowth, vp::ValueCoeffs;
                  capT::Int=10_000, nburn::Int=200, seed=42)
    Random.seed!(seed)  # Set specific seed
    shocks = randn(capT + nburn)

    return simulate(ncgm, vp, shocks; capT=capT, nburn=nburn)
end

"""
This function evaluates the Euler Equation residual for a single point (k, z)
"""
function EulerEquation!(out::Vector{Float64}, ncgm::NeoclassicalGrowth,
                        vp::ValueCoeffs, k::Float64, z::Float64,
                        nodes::Vector{Float64}, weights::Vector{Float64})
    # Evaluate consumption today
    k1 = env_condition_kp!(out, ncgm, vp, k, z)
    c = expendables_t(ncgm, k, z) - k1
    LHS = du(ncgm, c)

    # For each of realizations tomorrow, evaluate expectation on RHS
    RHS = 0.0
    for (eps, w) in zip(nodes, weights)
        # Compute ztp1
        z1 = z^ncgm.ρ * exp(eps)

        # Evaluate the ktp2
        ktp2 = env_condition_kp!(out, ncgm, vp, k1, z1)

        # Get c1
        c1 = expendables_t(ncgm, k1, z1) - ktp2

        # Update RHS of equation
        RHS = RHS + w*du(ncgm, c1)*(1 - ncgm.δ + df(ncgm, k1, z1))
    end

    return abs(ncgm.β*RHS/LHS - 1.0)
end

"""
Given simulations for k and z, it computes the euler equation residuals
along the entire simulation. It reports the mean and max values in
log10.
"""
function ee_residuals(ncgm::NeoclassicalGrowth, vp::ValueCoeffs,
                      ksim::Vector{Float64}, zsim::Vector{Float64}; Qn::Int=10)
    # Figure out how many periods we simulated for and make sure k and z
    # are same length
    capT = length(ksim)
    @assert length(zsim) == capT

    # Finer integration nodes
    eps_nodes, weight_nodes = qnwnorm(Qn, 0.0, ncgm.σ^2)
    temp = Array{Float64}(undef, n_complete(2, vp.d))

    # Compute EE for each period
    EE_resid = Array{Float64}(undef, capT)
    for t=1:capT
        # Pull out current state
        k, z = ksim[t], zsim[t]

        # Compute residual of Euler Equation
        EE_resid[t] = EulerEquation!(temp, ncgm, vp, k, z, eps_nodes, weight_nodes)
    end

    return EE_resid
end

function ee_residuals(ncgm::NeoclassicalGrowth, vp::ValueCoeffs; Qn::Int=10)
    # Simulate and then call other ee_residuals method
    ksim, zsim = simulate(ncgm, vp)

    return ee_residuals(ncgm, vp, ksim, zsim; Qn=Qn)
end

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractQuasiArray,B<:AbstractQuasiArray} =
    quasildivapplystyle(MemoryLayout(A), MemoryLayout(B))
ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractQuasiArray,B<:AbstractArray} =
    quasildivapplystyle(MemoryLayout(A), MemoryLayout(B))
ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractArray,B<:AbstractQuasiArray} =
    quasildivapplystyle(MemoryLayout(A), MemoryLayout(B))

quasildivapplystyle(_, _) = LdivApplyStyle()

similar(L::Ldiv, ::Type{T}, axes::Tuple{<:AbstractQuasiVector,<:AbstractQuasiVector}) where T = similar(QuasiArray{T}, axes)
similar(L::Ldiv, ::Type{T}, axes::Tuple{<:AbstractQuasiVector,<:AbstractVector}) where T = similar(QuasiArray{T}, axes)
similar(L::Ldiv, ::Type{T}, axes::Tuple{<:AbstractVector,<:AbstractQuasiVector}) where T = similar(QuasiArray{T}, axes)
similar(L::Ldiv, ::Type{T}, axes::Tuple{<:AbstractQuasiVector}) where T = similar(QuasiArray{T}, axes)

function copyto!(dest::QuasiArray, L::Ldiv{QuasiArrayLayout,QuasiArrayLayout})
    copyto!(parent(dest), Ldiv(Array(L.A), Array(L.B)))
    dest
end
function copyto!(dest::QuasiArray, L::Ldiv{QuasiArrayLayout})
    copyto!(parent(dest), Ldiv(Array(L.A), L.B))
    dest
end
function copyto!(dest::QuasiArray, L::Ldiv{<:Any,QuasiArrayLayout})
    copyto!(parent(dest), Ldiv(L.A, Array(L.B)))
    dest
end
copyto!(dest::AbstractArray, L::Ldiv{QuasiArrayLayout,QuasiArrayLayout}) =
    copyto!(dest, Ldiv(Array(L.A), Array(L.B)))
copyto!(dest::AbstractArray, L::Ldiv{QuasiArrayLayout}) =
    copyto!(dest, Ldiv(Array(L.A), L.B))
copyto!(dest::AbstractArray, L::Ldiv{<:Any,QuasiArrayLayout}) =
    copyto!(dest, Ldiv(L.A, Array(L.B)))

for op in (:pinv, :inv)
    @eval ApplyStyle(::typeof($op), args::Type{<:AbstractQuasiArray}) = QuasiArrayApplyStyle()
end
## PInvQuasiMatrix

const PInvQuasiMatrix{T, ARGS} = ApplyQuasiMatrix{T,typeof(pinv),ARGS}
const InvQuasiMatrix{T, ARGS} = ApplyQuasiMatrix{T,typeof(inv),ARGS}

parent(A::PInvQuasiMatrix) = first(A.args)
parent(A::InvQuasiMatrix) = first(A.args)
axes(A::PInvQuasiMatrix) = axes(Applied(A))
size(A::PInvQuasiMatrix) = map(length, axes(A))
pinv(A::PInvQuasiMatrix) = first(A.args)

@propagate_inbounds getindex(A::PInvQuasiMatrix{T}, k::Int, j::Int) where T =
    (Applied(A)*[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k]

*(A::PInvQuasiMatrix, B::AbstractQuasiMatrix, C...) = apply(*,Applied(A), B, C...)
*(A::PInvQuasiMatrix, B::MulQuasiArray, C...) = apply(*,Applied(A), Applied(B), C...)

## QuasiArray special case
inv(A::QuasiMatrix) = QuasiArray(inv(A.parent), reverse(A.axes))


text:
import Base: close, unsafe_read, unsafe_write

mutable struct FixedLengthRecord{T,C} <: Record
   io       :: IO    # underlying I/O stream
   reclen   :: T     # length of this record
   nleft    :: T     # bytes left in this record
   writable :: Bool  # whether this record is writable
   convert  :: C     # convert method
end

function Record( f::FortranFile{DirectAccess}, recnum::Integer )
## constructor for readable records
   conv = f.convert
   reclen = f.acctyp.reclen
   gotorecord(f, recnum)
   FixedLengthRecord(f.io, reclen, reclen, false, conv)
end

function Record( f::FortranFile{DirectAccess}, recnum::Integer, towrite::Integer )
## constructor for writable records
   conv = f.convert
   reclen = f.acctyp.reclen
   if towrite > reclen
      fthrow("attempting to write record of $(towrite) bytes into a file of record length $(reclen) bytes")
   end
   gotorecord(f, recnum)
   FixedLengthRecord(f.io, reclen, reclen, true, conv)
end

function gotorecord( f::FortranFile{DirectAccess}, recnum::Integer )
   reclen = f.acctyp.reclen
   pos = (recnum-1) * reclen
   seek(f.io, pos)
end

function unsafe_read( rec::FixedLengthRecord, p::Ptr{UInt8}, n::UInt )
   if (n > rec.nleft)
      fthrow("attempting to read beyond record end")
   end
   unsafe_read( rec.io, p, n )
   rec.nleft -= n
   nothing
end

function unsafe_write( rec::FixedLengthRecord, p::Ptr{UInt8}, n::UInt )
   if (n > rec.nleft); fthrow("attempting to write beyond record end"); end
   nwritten = unsafe_write( rec.io, p, n )
   rec.nleft -= n
   return nwritten
end

function close( rec::FixedLengthRecord )
   if rec.nleft != 0
      if rec.writable
         # Fortran standard 9.6.4.5.2 point 7:
         # "If the file is connected for direct access and the values specified by the
         #  output list do not fill the record, the remainder of the record is undefined."
         # Following gfortran, we fill it with zeros.
         write(rec.io, zeros(UInt8, rec.nleft))
      else
         skip(rec.io, rec.nleft)
      end
   end
   nothing
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using SymPy

function diffeq_3()
    @vars z t
    P, ∇Pk, Uk, dUdtk = symbols("P ∇Pk Uk dUdtk", cls = sympy.Function)
    P = cos(t)*cos(z)
    ∇Pk = diff.(P, z)
    dUdtk = diff.(Uk(z, t), t)
    diffeq = Eq(dUdtk, ∇Pk)
    display(diffeq)
    solution = pdsolve(diffeq)
end

diffeq_3()

# %%
using SymPy

function diffeq_3()
    @vars z t
    P, ∇Pk, Uk, dUdtk = symbols("P ∇Pk Uk dUdtk", cls = sympy.Function)
    P = cos(t)*cos(z)
    ∇Pk = diff.(P, z)
    dUdtk = diff.(Uk(z, t), t)
    diffeq = Eq(dUdtk, ∇Pk)
    solution = pdsolve(diffeq)
    [diffeq, solution]
end

diffeq_3()

# %%
@show diffeq_3();

# %%

max_stars_repo_path: 0016/SymPy example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/SymPy example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/SymPy example.jl

max_forks_repo_name: genkuroki/public


text:
module SymbolicSAT

using SymbolicUtils
using SymbolicUtils: Sym, Term, operation, arguments, Symbolic, symtype, istree
using Z3

export Constraints, issatisfiable, isprovable

# Symutils -> Z3

function to_z3(term, ctx)
    if istree(term)
        op = operation(term)
        args = arguments(term)

        # weird special case
        if length(args) == 1 && (op === (!) || op === (~))
            op = not
        elseif op === (&)
            op = and
        elseif op === (|)
            op = or
        end
        args′ = map(x->to_z3(x, ctx), args)
        s = findfirst(x->x isa Symbolic, args′)
        if s !== nothing
            error("$s of symtype $(symtype(s)) was not converted into a z3 expression")
        end
        op(args′...)
    else
        term
    end
end

for (jlt, z3t) in [Integer => :int, Real => :real]
    @eval to_z3(x::Sym{<:$jlt}, ctx) = $(Symbol(z3t, "_const"))(ctx, String(nameof(x)))
end

struct Constraints
    constraints::Vector
    solver::Z3.Solver
    context::Z3.Context

    function Constraints(cs::Vector, solvertype="QF_NRA")
        ctx = Context()
        s = Solver(ctx, solvertype)

        for c in cs
            add(s, to_z3(c, ctx))
        end

        new(cs, s, ctx)
    end
end

function Base.show(io::IO, c::Constraints)
    cs = c.constraints
    println("Constraints:")
    for i in 1:length(cs)
        print(io, "  ")
        print(io, cs[i])
        i != length(cs) && println(io, " ∧")
    end
end

function issatisfiable(expr::Symbolic{Bool}, cs::Constraints)
    Z3.push(cs.solver)
    add(cs.solver, to_z3(expr, cs.context))
    res = check(cs.solver)
    Z3.pop(cs.solver,1)
    if res == Z3.sat
        return true
    elseif res == Z3.unsat
        return false
    elseif res == Z3.unknown
        return nothing
    end
end

function isprovable(expr, cs::Constraints)
    sat = issatisfiable(expr, cs)
    sat === true ? !issatisfiable(!expr, cs) : false
end

issatisfiable(expr::Bool, Constraints) = expr

boolsym(x::Symbolic) = symtype(x) == Bool
boolsym(x) = false

function resolve(x, ctx)
     isprovable(x, ctx) === true ?
        true : isprovable(!(x), ctx) === true ? false : x
end

end # module

max_stars_repo_name: JuliaSymbolics/SymbolicSAT

max_issues_repo_name: JuliaSymbolics/SymbolicSAT

max_forks_repo_name: JuliaSymbolics/SymbolicSAT


text:
# avoid some images failing on travis
_backend_skips = Dict{Symbol, Vector{Int}}(
    :gr => [],
    :pyplot => [],
    :plotlyjs => [24],
    :pgfplots => [],
)

function no_initialize_backend(be)
    Plots.CURRENT_BACKEND.sym = be
    Plots.CURRENT_BACKEND.pkg = Plots._backend_instance(be)
end

# ------------------------------------------------------------------------------
# reference images

filename(i::Int) = string("ref", i, i in Plots._animation_examples ? ".gif" : ".png")

"""
    generate_reference_images(be::Symbol, overwrite = true)

Generate all Plots.jl reference images for the backend `be` for PlotDocs.jl.
Working backends are `:gr`, `:pyplot`, `:plotlyjs` and `:pgfplots`.
If `overwrite == true` existing files are overwritten.
"""
function generate_reference_images(be::Symbol, overwrite = true)
    for i in eachindex(Plots._examples)
        if i in Plots._backend_skips[be] || i in PlotDocs._backend_skips[be]
            @info "Skipping $be reference image $i."
        else
            @info "Generating $be reference image $i."
            generate_reference_image(be, i, overwrite)
        end
    end
end

"""
    function generate_reference_image(be::Symbol, i::Int, overwrite::Bool = true)

Generate the Plots.jl reference image `i` for the backend `be` for PlotDocs.jl.
Working backends are `:gr`, `:pyplot`, `:plotlyjs` and `:pgfplots`.
If `overwrite == true` the existing file is overwritten.
"""
function generate_reference_image(be::Symbol, i::Int, overwrite = true)
    @eval PlotDocs theme(:default)
    no_initialize_backend(be)
    map(ex -> Base.eval(PlotDocs, ex), Plots._examples[i].exprs)

    dir = joinpath(IMGDIR, string(be))
    isdir(dir) || mkpath(dir)

    fn = joinpath(dir, filename(i))
    overwrite || !isfile(fn) || return

    if i in REFERENCE_ANIMATION_INDS
        anim = @eval PlotDocs anim
        Plots.gif(anim, fn, fps = 15)
    else
        Plots.savefig(fn)
    end

    return
end

const REFERENCE_ANIMATION_INDS = (2, 30)

# ------------------------------------------------------------------------------
# documentation examples

"""
    generate_doc_images(overwrite = true)

Generate all PlotDocs images.
If `overwrite == true` existing files are overwritten.
"""
function generate_doc_images(overwrite = true)
    for (i, id) in enumerate(keys(DOC_IMAGES))
        @info "Generating PlotDocs image $i: $id"
        generate_doc_image(id, overwrite)
    end
end

"""
    function generate_doc_image(id::AbstractString, overwrite::Bool = true)

Generate the PlotDocs image specified by `id`.
Choose from $(keys(DOC_IMAGES)).
If `overwrite == true` the existing file is overwritten.
"""
function generate_doc_image(id::AbstractString, overwrite = true)
    @eval PlotDocs theme(:default)
    Random.seed!(1234)
    map(ex -> Base.eval(PlotDocs, ex), DOC_IMAGES[id])

    dir = joinpath(IMGDIR, DOC_IMAGE_FILES[id])
    isdir(dir) || mkpath(dir)

    filetype = if id in DOC_ANIMATION_IMAGES
        ".gif"
    elseif id in DOC_MP4_IMAGES
        ".mp4"
    else
        ".png"
    end
    fn = joinpath(dir, string(id, filetype))
    overwrite || !isfile(fn) || return

    if id in DOC_ANIMATION_IMAGES
        anim = @eval PlotDocs anim
        if typeof(anim) <: Plots.AnimatedGif
            mv(anim.filename, fn, force = true)
        else
            Plots.gif(anim, fn, fps = 15)
        end
    elseif id in DOC_MP4_IMAGES
        anim = @eval PlotDocs anim
        mp4(anim, fn, fps = 30)
    else
        Plots.savefig(fn)
    end

    return
end


text:
#*************************************************************************#
#
# Abstract Types
#
#*************************************************************************#
"""
    ParametricModel

Define a static model relating input and output.
"""
abstract type StaticModel end

"""
    ParametricModel

Define a static parametric model `z = f(x, Θ)`.

Basic Functions
---------------

| Function                          | Description                                                 |
|:--------------------------------- |:----------------------------------------------------------- |
| `evaluate!(f, x, Θ, z[, dx, dΘ])` | Compute the output `z` and its first derivatives            |
| `z_buffer(f)`                     | Generate `z` buffer to be used in `evaluate!`               |
| `dx_buffer(f)`                    | Generate `dx` buffer to be used in `evaluate!`              |
| `dΘ_buffer(f)`                    | Generate `dΘ` buffer to be used in `evaluate!`              |
| `isparametric(f)`                 | Return `true` if the given StaticModel is a ParametricModel |
| `length_input(f)`                 | Return number of model inputs                               |
| `length_output(f)`                | Return number of model outputs                              |
| `length_params(f)`                | Return number of model parameters                           |
"""
abstract type ParametricModel <: StaticModel end

"""
    AparametricModel

Define a static aparametric model `z = f(x)`.

Basic Functions
---------------

| Function                   | Description                                                   |
|:-------------------------- |:-----------------------------------------------------------   |
| `evaluate!(f, x, z[, dx])` | Compute the output `z` and its first derivatives              |
| `z_buffer(f)`              | Generate `z` buffer to be used in `evaluate!`                 |
| `dx_buffer(f)`             | Generate `dx` buffer to be used in `evaluate!`                |
| `isparametric(f)`          | Return `false` if the given StaticModel is a AparametricModel |
| `length_input(f)`          | Return number of model inputs                                 |
| `length_output(f)`         | Return number of model outputs                                |
"""
abstract type AparametricModel <: StaticModel end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
"""
    z_buffer(f::StaticModel) -> z

Return the array to be used as buffer for the output vector `z`
"""
function z_buffer(f::StaticModel)
    return Vector{Float64}(f.noutputs)
end

"""
    dx_buffer(f::StaticModel) -> dx

Return the array to be used as buffer for the Jacobian matrix `dx`.
"""
function dx_buffer(f::StaticModel)
    return Matrix{Float64}(f.noutputs, f.ninputs)
end

"""
    dΘ_buffer(f::ParametricModel) -> dΘ

Return the array to be used as buffer for the Jacobian matrix `dΘ`.
"""
function dΘ_buffer(f::ParametricModel)
    return Matrix{Float64}(f.noutputs, f.nparams)
end

#*************************************************************************#
#
# Compute Output and its Derivatives
#
#*************************************************************************#
"""
    evaluate!(f::ParametricModel, x, Θ, z[, dx, dΘ])

For a static function `z = f(x, Θ)` compute the
output `z` and its first derivatives `dx` and `dΘ.

The real vectors `x` and `Θ` contains respectively the
input and parameters to the static model `f`.  The computed
values are written in `z`, `dx` and `dΘ`. Which are the static
function output and the jacobian matrices containing derivatives in
relation to `x` and `Θ`.
"""
function evaluate!(f::ParametricModel,
                   x::AbstractVector{<:Real},
                   Θ::AbstractVector{<:Real},
                   z::AbstractVector{Float64};
                   dx::AbstractMatrix{Float64}=Matrix{Float64}(0, 0),
                   dΘ::AbstractMatrix{Float64}=Matrix{Float64}(0, 0))

    # ---
    # Check Arguments
    # ---
    nx = length(x)
    nparams = length(Θ)
    nz = length(z)
    size_dx = size(dx)
    size_dΘ = size(dΘ)

    if nx != f.ninputs
        throw(ArgumentError(string("Vector `x` have $nx elements when it should ",
                                   "have $(f.ninputs).")))
    end
    if nparams != f.nparams
        throw(ArgumentError(string("Vector `Θ` have $nparams elements when it ",
                                   "should have $(f.nparams).")))
    end
    if nz != f.noutputs
        throw(ArgumentError(string("Vector `z` have $nz elements when it should ",
                                   "have $(f.noutputs).")))
    end
    if !isempty(dx) && size_dx != (f.noutputs, f.ninputs)
        throw(ArgumentError(string("The matrix `dx` have dimension $(size_dx) ",
                                   "when it should have dimension ",
                                   "$((f.noutputs, f.ninputs)).")))
    end
    if !isempty(dΘ) && size_dΘ != (f.noutputs, f.nparams)
        throw(ArgumentError(string("The matrix `dΘ` have ",
                                   "dimension $(size_dΘ) ",
                                   "when it should have dimension ",
                                   "$((f.noutputs, f.nparams)).")))
    end
    # ---

    # Evaluate
    _evaluate!(f, x, Θ, z, dx, dΘ)
end

"""
    evaluate!(f::ParametricModel, x, z[, dx])

For a static function `z = f(x)` compute the
output `z` and its first derivatives `dx`.

The real vector `x` contain respectively the
input and parameters to the static model `f`.  The computed
values are written in `z`, `dx`. Which are the static
function output and the  jacobian matrix.
"""
function evaluate!(f::AparametricModel,
                   x::AbstractVector{<:Real},
                   z::AbstractVector{Float64};
                   dx::AbstractMatrix{Float64}=Matrix{Float64}(0, 0))

    # ---
    # Check Arguments
    # ---
    nx = length(x)
    nz = length(z)
    size_dx = size(dx)

    if nx != f.ninputs
        throw(ArgumentError(string("Vector `x` have $nx elements when it should ",
                                   "have $(f.ninputs).")))
    end
    if nz != f.noutputs
        throw(ArgumentError(string("Vector `z` have $nz elements when it should ",
                                   "have $(f.noutputs).")))
    end
    if !isempty(dx) && size_dx != (f.noutputs, f.ninputs)
        throw(ArgumentError(string("The matrix `dx` have dimension $(size_dx) ",
                                   "when it should have dimension ",
                                   "$((f.noutputs, f.ninputs)).")))
    end
    # ---

    # Evaluate
    _evaluate!(f, x, z, dx)
end

#*************************************************************************#
#
# Auxiliar Functions
#
#*************************************************************************#
function isparametric(f::StaticModel)
    typeof(f) <: ParametricModel
end

function length_input(f::StaticModel)
    return f.ninputs
end

function length_output(f::StaticModel)
    return f.noutputs
end

function length_params(f::ParametricModel)
    return f.nparams
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module SubModule

buildfile = joinpath(@__DIR__, "..", "deps", "buildartifact")
if !isfile(buildfile)
    error("Package built incorrectly")
else
    include(buildfile)
end

f() = 1

end # module

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card(
    [
        dbc_cardheader("This is the header"),
        dbc_cardbody([
            html_h4("Card title", className="card-title"),
            html_p("This is some card text", className="card-text"),
        ]),
        dbc_cardfooter("This is the footer"),
    ],
    style=Dict("width" => "18rem"),
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module juliaAllotter

using CSV
using LinearAlgebra
using Statistics
using JLD
include("core.jl")

import Random
Random.seed!(2^13 - 1)

courseFile = "data/courseFile.csv"
studentsFile = "data/studentsFile.csv"
choiceIndices = [0, 1, 2, 3, 4, 5]
choiceWeights = [10, 1, 2, 3, 4, 5]
costWeights = [-1.0, 1.0, 1.0]
iterations = 4900

coursesDF = CSV.read(courseFile)
studentsDF = CSV.read(studentsFile)

courseCount = length(coursesDF[:, 1])
studentCount = length(studentsDF[:, 1])
cpiArray = Array(studentsDF.CPI)
times = Array(coursesDF.CourseNeeds)

choiceIdx = makeArray(studentCount, studentsDF, courseCount, coursesDF, choiceIndices)
choiceWeights = makeArray(studentCount, studentsDF, courseCount, coursesDF, choiceWeights)

initialAllotment = allotment(studentCount, courseCount, times, Matrix{Float64}(I, studentCount, studentCount))
initialAllotment = squeezeAllotment(initialAllotment)

choiceGoodnessOld, cpiGoodnessOld = calcGoodness(initialAllotment, choiceWeights, cpiArray)

finalAllottment, utility = runMCMC(initialAllotment, iterations, studentCount, courseCount, costWeights, choiceIdx, choiceWeights, studentsDF, cpiArray, times)
choiceGoodnessNew, cpiGoodnessNew = calcGoodness(finalAllottment, choiceWeights, cpiArray)

writePerformance(finalAllottment, choiceGoodnessOld, cpiGoodnessOld, choiceGoodnessNew, cpiGoodnessNew, utility)

end # module

max_stars_repo_name: saurabhkm/juliaAllotter

max_issues_repo_name: saurabhkm/juliaAllotter

max_forks_repo_name: saurabhkm/juliaAllotter


text:
struct MultiFieldCellField{R,S} <: CellField
  single_fields::Vector{<:CellField}
  cell_map::AbstractArray
  ref_trait::Val{R}
  cell_axes::AbstractArray
  metasize::Val{S}
  function MultiFieldCellField(single_fields::Vector{<:CellField})
    @assert length(single_fields) > 0
    f1 = first(single_fields)
    cell_map = get_cell_map(f1)
    ref_style = RefStyle(f1)
    cell_axes = get_cell_axes(f1)
    metasize = MetaSizeStyle(f1)
    R = get_val_parameter(ref_style)
    S = get_val_parameter(metasize)
    new{R,S}(single_fields,cell_map,ref_style,cell_axes,metasize)
  end
end

Arrays.get_array(a::MultiFieldCellField) = @notimplemented

CellData.get_memo(a::MultiFieldCellField) = @notimplemented

CellData.get_cell_map(a::MultiFieldCellField) = a.cell_map

CellData.get_cell_axes(a::MultiFieldCellField) = a.cell_axes

CellData.RefStyle(::Type{<:MultiFieldCellField{R}}) where R = Val(R)

CellData.MetaSizeStyle(::Type{<:MultiFieldCellField{R,S}}) where {R,S} = Val(S)

Base.length(f::MultiFieldCellField) = length(first(f.single_fields))

num_fields(a::MultiFieldCellField) = length(a.single_fields)

Base.getindex(a::MultiFieldCellField,i) = a.single_fields[i]

Base.iterate(a::MultiFieldCellField)  = iterate(a.single_fields)

Base.iterate(a::MultiFieldCellField,state)  = iterate(a.single_fields,state)

function Fields.evaluate(cf::MultiFieldCellField,x::AbstractArray)
  s = """
  Evaluation of MultiFieldCellField objects is not implemented at this moment.
  You need to extract the individual fields and then evaluate them separatelly.
  If ever implemented, evaluating a `MultiFieldCellField` directly would provide,
  at each evaluation point, a tuple with the value of the different fields.
  """
  @notimplemented s
end

function Geometry.restrict(a::MultiFieldCellField,trian::Triangulation)
  f = (ai) -> restrict(ai,trian)
  blocks = map(f,a.single_fields)
  blocks
end


text:
module NewtonsMethod

using ForwardDiff

function newtonroot(f, fp; x0, tol=1e-7, maxiter=1000)
    abserror = Inf
    iter = 1
    x = x0
    while abserror > tol && iter <= maxiter
        x_new = x - f(x)/fp(x)
        iter = iter +1
        abserror = abs(x_new - x)
        x = x_new
    end
    return x
end

function newtonroot(f; x0, tol = 1e-7, maxiter=1000)
    fp = x -> ForwardDiff.derivative(f, x)
    newtonroot(f, fp; x0=x0, tol=tol, maxiter=maxiter)
end

export newtonroot
end


text:
using BinaryBuilder

name = "LibUV"
version = v"2"

# Collection of sources required to build libuv
sources = [
    GitSource("https://github.com/JuliaLang/libuv.git",
              "4726f5ef94fe5e11335a6d6ef048b00a72110135"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libuv/

# Touch some files so that the build system doesn't try to re-run `acreconf`:
touch -c aclocal.m4
touch -c Makefile.in
touch -c configure

# `--with-pic` isn't enough; we really really need -fPIC and -DPIC everywhere...
# everywhere, especially on FreeBSD
./configure --prefix=$prefix --host=$target --with-pic CFLAGS="${CFLAGS} -DPIC -fPIC" CXXFLAGS="${CXXFLAGS} -DPIC -fPIC"
make -j${nproc} V=1
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libuv", :libuv)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# This file is part of Intersections.jl
export imerge2

"""
    imerge2(A, B, output=eltype(A)[])

Computes the intersection of `A` and `B` (sorted arrays), using a merge-like algorithm, and stores it in `output`.
"""
function imerge2(A, B, output=eltype(A)[])
    i = j = 1
    n, m = length(A), length(B)

    @inbounds while i <= n && j <= m
        c = cmp(_get_key(A, i), _get_key(B, j))
        if c == 0
            push!(output, _get_key(A, i))
            i += 1
            j += 1
        elseif c < 0
            i += 1
        else
            j += 1
        end
    end

    output
end


text:
module ImageComponentAnalysis

using Images, IndexedTables, JuliaDBMeta, OffsetArrays, StaticArrays, DataStructures

# Used in generic_labelling.jl to allow @nexprs macros.
using Base.Cartesian

abstract type MeasurementProperties end

struct Measurement <: MeasurementProperties
	area::Bool
	perimeter::Bool
end
Measurement(; area::Bool = true, perimeter::Bool = true) = Measurement(area, perimeter)

struct BoundingBox <: MeasurementProperties
	box_area::Bool
end
BoundingBox(; box_area::Bool = true) = BoundingBox(box_area)

struct BasicTopology <: MeasurementProperties
	holes::Bool
	euler_number::Bool
end
BasicTopology(; holes::Bool = true, euler_number::Bool = true) = BasicTopology(holes, euler_number)

struct RegionEllipse <: MeasurementProperties
	centroid::Bool
	semi_axes::Bool
	orientation::Bool
	eccentricity::Bool
end
RegionEllipse(; semi_axes::Bool = true, orientation::Bool = true, eccentricity::Bool = true, centroid::Bool = true) = RegionEllipse(centroid, semi_axes, orientation, eccentricity)

struct Boundary <: MeasurementProperties
	inner::Bool
	outer::Bool
	#holes_inner::Bool
	#holes_outer::Bool
end
Boundary(; inner::Bool = true, outer::Bool = true) = Boundary(inner, outer)

abstract type ComponentLabellingAlgorithm end
struct OneComponent2D <: ComponentLabellingAlgorithm end
struct OneComponent3D <: ComponentLabellingAlgorithm end
struct Generic <: ComponentLabellingAlgorithm end
struct ContourTracing <: ComponentLabellingAlgorithm end

abstract type TracingAlgorithm end
struct CostaOuter <: TracingAlgorithm end
struct MooreInner <: TracingAlgorithm end

abstract type Connectivity end
struct FourConnected <: Connectivity end
struct EightConnected <: Connectivity end

include("boundingbox.jl")
include("measurement.jl")
include("basic_topology.jl")
include("region_ellipse.jl")
include("boundary.jl")
include("measure_components.jl")
include("contour_tracing.jl")
include("costa_outer.jl")
include("moore_inner.jl")
include("one_component_2d.jl")
include("one_component_3d.jl")
include("generic_labelling.jl")

export
	Boundary,
	BoundingBox,
	Measurement,
	BasicTopology,
	RegionEllipse,
	measure_components,
	label_components,
	ContourTracing,
	FourConnected,
	EightConnected,
	OneComponent2D,
  OneComponent3D,
	Generic,
	trace_boundary,
	CostaOuter,
	MooreInner
end # module


text:
function fermat_test(n::Int128)
    for k in 1:100
        a::Int128 = rand(2:n-1)

        if gcd(n, a) != 1
            return false
        end
        if a^(n-1) % n != 1
            return false
        end
    end

    true
end

println("57 is prime: ", fermat_test(Int128(57)))
println("23 is prime: ", fermat_test(Int128(23)))


text:
"""
---------------------------------------------------
Types and Constructors
---------------------------------------------------
"""

mutable struct VERTICAL_DESCENT_SAMPLER <: SAMPLER
	z_init::Sampleable
	z_descent::Sampleable
	vx_init::Sampleable
	vz_init::Sampleable
	θ::Sampleable
	t_extra::Sampleable
	t_descent::Sampleable
end

function vertical_descent_sampler(;z_init = Uniform(300,600),
						z_descent = Normal(100,8),
						vx_init = Uniform(30,45),
						vz_init = Uniform(-500fpm2fps,-300fpm2fps),
						θ = Truncated(Normal(6,3),2,15),
						t_extra = Uniform(0.1,0.3), #Normal(10,3),
						t_descent = Normal(30,5))
	return VERTICAL_DESCENT_SAMPLER(z_init, z_descent, vx_init, vz_init, θ, t_extra, t_descent)
end

mutable struct VERTICAL_DESCENT <: UAM_TRAJECTORY
	z_init::Float64
	z_descent::Float64
	vx_init::Float64
	vz_init::Float64
	θ::Float64
	t_extra::Float64
	t_descent::Float64
	sampler::SAMPLER
	dt::Float64
	vz_max::Float64
	az_max::Float64
	ax_max::Float64
	p::Array{Float64,2}
	v::Array{Float64,2}
end

function vertical_descent(;z_init = 0.0,
						z_descent = 0.0,
						vx_init = 0.0,
						vz_init = 0.0,
						θ = 0.0,
						t_extra = 0.0,
						t_descent = 0.0,
						sampler = vertical_descent_sampler(),
						dt = 1.0,
						vz_max = 550fpm2fps,
						az_max = 0.3g_ft,
						ax_max = 0.1g,
						p = zeros(1,3),
						v = zeros(1,3))
	return VERTICAL_DESCENT(z_init, z_descent, vx_init, vz_init, θ, t_extra, t_descent, sampler, dt,
							 vz_max, az_max, ax_max, p, v)
end

"""
---------------------------------------------------
Functions
---------------------------------------------------
"""

function solve_trajectory!(τ::VERTICAL_DESCENT)
	# Glide portion #################################################
	# Get some needed info
	τ.vz_init = clamp(τ.vz_init, -τ.vz_max, Inf)
	ttot = abs(τ.z_init/τ.vz_init) + τ.t_extra*abs(τ.z_init/τ.vz_init)
	N = convert(Int64, floor(ttot/τ.dt))
	# Difference matrix
	D = zeros(N-1,N)
	for i = 1:(N-1)
	    D[i,i] = -1
	    D[i,i+1] = 1
	end
	D = D./τ.dt

	# Create and solve the optimization problem
	px = Variable(N)
	pz = Variable(N)
	vx = Variable(N)
	vz = Variable(N)

	objec = norm(D*vx) + λ*norm(vz*ft2m + vx*tand(τ.θ))

	constraints = [[px[i+1] == px[i] + vx[i]*τ.dt for i = 1:N-1];
				   [pz[i+1] == pz[i] + vz[i]*τ.dt for i = 1:N-1];
				   vx[1] == τ.vx_init; px[1] == 0;
				   vz[1] == τ.vz_init; pz[1] == τ.z_init;
				   vx[N] == 0; pz[N] == τ.z_descent; vz[N] == 0;
				   vx ≥ 0; pz ≥ 0; vz ≤ 0; vz ≥ -τ.vz_max;
				   D*vx ≤ τ.ax_max; D*vz*ft2m ≤ τ.az_max]

	prob1 = minimize(objec, constraints)
	solve!(prob1, ECOS.Optimizer(verbose=0))

	# Vertical descent portion ######################################
	Nd = convert(Int64, floor(τ.t_descent/τ.dt))
	# Difference matrix
	Dd = zeros(Nd-1,Nd)
	for i = 1:(Nd-1)
	    Dd[i,i] = -1
	    Dd[i,i+1] = 1
	end
	Dd = Dd./τ.dt

	pzd = Variable(Nd)
	vzd = Variable(Nd)

	objec = norm(Dd*vzd)

	constraints = [[pzd[i+1] == pzd[i] + vzd[i]*τ.dt for i = 1:Nd-1];
				   vzd[1] == 0; pzd[1] == τ.z_descent;
				   pzd[Nd] == 0; vzd[Nd] == 0;
				   pzd ≥ 0; vzd ≤ 0; vzd ≥ -τ.vz_max;
				   Dd*vzd*ft2m ≥ -τ.az_max; Dd*vzd*ft2m ≤ τ.az_max]

	prob2 = minimize(objec, constraints)
	solve!(prob2, ECOS.Optimizer(verbose=0))

	# Update the trajectory
	p = zeros(N+Nd,3)
	p[1:N,1] = px.value[:,1]
	p[N+1:end,1] = px.value[end,1]*ones(Nd)
	p[1:N,3] = pz.value[:,1]
	p[N+1:end,3] = pzd.value
	v = zeros(N+Nd,3)
	v[1:N,1] = vx.value[:,1]
	v[1:N,3] = vz.value[:,1]
	v[N+1:end,3] = vzd.value

	τ.p = copy(p)
	τ.v = copy(v)

	return Int(prob1.status) == 1
end

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


text:
using JuLIP, ACE, JuLIP.Potentials, LinearAlgebra
using ACE.Testing: lsq, get_V0
using LinearAlgebra: qr, cond
using Plots

using ACE
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d
using JuLIP.Potentials: i2z, numz
using ACE.RPI
using ACEds
using ACEds.equivRotations3D
using ACEds.equivRotations3D: _mrange1, ERot3DCoeffs
#using ACE.SphericalHarmonics: *
using ACE.SphericalHarmonics
using ACE.SphericalHarmonics: PseudoSpherical, index_y, cart2spher, sizeY, sizeP, allocate_p, compute_coefficients, compute_p!, cYlm!

using JuLIP.MLIPs

using ReferenceFrameRotations: angle_to_dcm

using StaticArrays

using Random

## Libraries added by Matthias

using Combinatorics
using StatsBase
##

function rand_unit_vec(d::Int64)
	u = [rand() for i in 1:d ]
	u/=norm(u)
	return SVector{d}(u)
end

N = 3 # N-body interaction
r_vec = [rand_unit_vec(3) for i in 1:N]

"""
Specify parameters
"""
ll = @SVector [2,1,3]
@assert N  == length(ll)
L = maximum(ll)

coeff = compute_coefficients(L) # Precompute coefficients for evaluation of spherical Harmonics
"""
Precompute/evaluate complex spherical harmonics
"""
Y = [ Array{ComplexF64}(undef,sizeY(L)) for i in 1:N]
let P =  Array{Float64}(undef, sizeP(L)) # used only temporarily
	for i in 1:N
		S = cart2spher(r_vec[i])
		compute_p!(L, S, coeff, P)
		cYlm!(Y[i], L, S, P)
	end
end

"""Define rotation"""
Q = angle_to_dcm(.5,.4,-.3, :ZYX)

"""Precompute/evaluate rotated complex spherical harmonics"""
Y_rot = [ Array{ComplexF64}(undef,sizeY(L)) for i in 1:N]
let P =  Array{Float64}(undef, sizeP(L))
	for i in 1:N
		S = cart2spher(Q*r_vec[i])
		compute_p!(L, S, coeff, P)
		ACE.SphericalHarmonics.cYlm!(Y_rot[i], L, S, P)
	end
end

"""Precompute coefficients for rotational equivariant basis"""
erotc = ERot3DCoeffs(Float64)

function reB(Y, k::Int, ll::SVector, mm::SVector, A::ERot3DCoeffs{T}) where {T}
	"""
	compute symmetrized precomputed values Y of spherical harmonics
	"""
	N  = length(Y)
	@assert length(Y) == length(ll)
	L = maximum(ll)

	#Y_vec = [prod([ Y[j][index_y(ll[j],mu[j])] for j in 1:N ])
	#	for  (i,mu) in enumerate(_mrange1(ll))]

	eB_val = sum(
		[A(ll,mu,mm)[:,k] * prod([ Y[j][index_y(ll[j],mu[j])] for j in 1:N ])
		for  (i,mu) in enumerate(_mrange1(ll))]
	)

	return eB_val
end

function rpiB(Y, k::Int, ll::SVector, mm::SVector, A::ERot3DCoeffs{T}) where {T}
	"""
	compute symmetrized precomputed values Y of spherical harmonics
	"""
	n  = length(Y)
	@assert length(Y) == length(ll)
	L = maximum(ll)

	#Y_vec = [prod([ Y[j][index_y(ll[j],mu[j])] for j in 1:N ])
	#	for  (i,mu) in enumerate(_mrange1(ll))]
	erpiB_val = zeros(ComplexF64, 3)
	@assert length(Y) == length(ll)
	for p in Combinatorics.permutations(1:n)
		rpi_basis_val += reBreB(Y, k, ll, mm, A)
	end

	return erpiB_val/factorial(n)
end

#allocate_p(L)

#mm = @SVector [-1,1]
#eB_vec = [eB(Y, k, ll, mm, erotc ) for  (i,mm) in enumerate(_mrange1(ll))]
#eB_rot_vec = [eB(Y_rot, k, ll, mm, erotc ) for  (i,mm) in enumerate(_mrange1(ll))]
#print(eBk,"\n")
tol = 1E-10
for k in 1:3
	difference_vec = [Q*eB(Y, k, ll, mm, erotc ) - eB(Y_rot, k, ll, mm, erotc ) for  (i,mm) in enumerate(_mrange1(ll))]
	print("k = ", k, " ", all([norm(dd) < tol for dd in difference_vec]),"\n")
end
	#print([norm(dd) for dd in difference_vec])

for (im, (jm,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll))), (ik, (jk,kk)) in enumerate(Iterators.product(1:3, _mrange1(ll)))
	print(im,", ", ik,"\n")
end

dd  = 3
function re_gramian(A::ERot3DCoeffs{T}, ll::SVector) where {T}
	len = length(_mrange1(ll))*3
	GG = zeros(ComplexF64, len, len) # = Gramian of rotational eqiuvariant span
	# kk = mprime
    for (i, (k,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll))), (ip, (kp,mmp)) in enumerate(Iterators.product(1:3, _mrange1(ll)))
		for mu in _mrange1(ll)
		  	GG[i,ip] += dot(A(ll, mm, mu)[:,k],A(ll, mmp, mu)[:,kp])
		end
	end
	return GG
end

GG = re_gramian(erotc, ll)

"""return coefficients for a rotational invariant tensor basis"""
function re_basis(A::ERot3DCoeffs, ll::SVector{N, Int}) where {N}
	#Mre = collect(Iterators.product(1:3, _mrange1(ll)))
	GG = re_gramian(A, ll)
   	S = svd(GG)
   	rk = rank(GG; rtol =  1e-7)
	Ure = S.U[:, 1:rk]'

	return Diagonal(1. ./ sqrt.(S.S[1:rk])) * Ure
end

re_coeff  = re_basis(erotc, ll)
print("Basis size: ", size(re_coeff))

function eval_re_basis(re_coeff,Y, A::ERot3DCoeffs, ll::SVector{N, Int}) where {N}
	b_size, span_size  = size(re_coeff)
	re_basis_val = zeros(ComplexF64, 3, b_size)
	eB_vec = zeros(ComplexF64, 3, span_size)
	for (ii, (k,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll)))
		eB_vec[:,ii] = eB(Y, k, ll, mm, A)
	end	#print(size(eB_vec))
	for i in 1:b_size
		re_basis_val[:,i] = sum([re_coeff[i,ii] * eB_vec[:,ii] for (ii, (k,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll)))])
	end
	return re_basis_val
end

function rpi_basis(A::Rot3DCoeffs,
						 zz::SVector{N},
						 nn::SVector{N, Int},
						 ll::SVector{N, Int}) where {N}
	Ure = re_basis(A, ll)
	Mre = collect( Iterators.product(1:3, _mrange1(ll)) )   # rows...
	G = _gramian(zz, nn, ll, Ure, Mre)
   	S = svd(G)
   	rk = rank(G; rtol =  1e-7)
	Urpi = S.U[:, 1:rk]'
	return Diagonal(sqrt.(S.S[1:rk])) * Urpi * Ure, Mre
end

function rpi_gramian_old(A::ERot3DCoeffs{T}, ll::SVector) where {T}
	len = length(_mrange1(ll))*3
	n = length(ll)
	GG = zeros(ComplexF64, len, len) # = Gramian of rotational eqiuvariant span
	# kk = mprime
    for (i, (k,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll))), (ip, (kp,mmp)) in enumerate(Iterators.product(1:3, _mrange1(ll)))
		for mu in _mrange1(ll)
			for perm in Combinatorics.permutations(1:n)
		  		GG[i,ip] += dot(A(ll[perm], mm[perm], mu[perm])[:,k],A(ll[perm], mmp[perm], mu[perm])[:,kp])
			end
		end
	end
	return GG
end

function rpi_gramian(zz, nn, ll, Ure, Mre)
   N = length(nn)
   nre = size(Ure, 1)
   @assert size(Ure, 1) == nre
   G = zeros(nre, nre)
   for σ in permutations(1:N)
      if (zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mre), (iU2, mm2) in enumerate(Mre)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Ure[i1, iU1]) * Ure[i2, iU2]
            end
         end
      end
   end
   return G
end

function ri_gramian(zz, nn, ll, Uri, Mri)
   N = length(nn)
   nri = size(Uri, 1)
   @assert size(Uri, 1) == nri
   G = zeros(nri, nri)
   for σ in permutations(1:N)
      if (zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mri), (iU2, mm2) in enumerate(Mri)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Uri[i1, iU1]) * Uri[i2, iU2]
            end
         end
      end
   end
   return G
end

rpi_basis(A::Rot3DCoeffs, zz, nn, ll) =
			rpi_basis(A, SVector(zz...), SVector(nn...), SVector(ll...))

function rpi_basis(A::Rot3DCoeffs,
						 zz::SVector{N},
						 nn::SVector{N, Int},
						 ll::SVector{N, Int}) where {N}
	Uri = ri_basis(A, ll)
	Mri = collect( _mrange(ll) )   # rows...
	G = _gramian(zz, nn, ll, Uri, Mri)
   	S = svd(G)
   	rk = rank(G; rtol =  1e-7)
	Urpi = S.U[:, 1:rk]'
	return Diagonal(sqrt.(S.S[1:rk])) * Urpi * Uri, Mri
end

function _gramian(zz, nn, ll, Uri, Mri)
   N = length(nn)
   nri = size(Uri, 1)
   @assert size(Uri, 1) == nri
   G = zeros(nri, nri)
   for σ in permutations(1:N)
      if (zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mri), (iU2, mm2) in enumerate(Mri)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Uri[i1, iU1]) * Uri[i2, iU2]
            end
         end
      end
   end
   return G
end

function eval_rpi_basis(re_coeff,Y, A::ERot3DCoeffs, ll::SVector{N, Int}) where {N}

	rpi_basis_val = zeros(ComplexF64, 3, size(re_coeff)...)
	n  = length(Y)
	@assert length(Y) == length(ll)
	for p in Combinatorics.permutations(1:n)
		rpi_basis_val += eval_re_basis(re_coeff,Y[p], A, ll)
	end
	re_basis_val/=factorial(n)
	return rpi_basis_val
end

eB_basis  = eval_re_basis(re_coeff, Y, erotc, ll)
eB_basis_rot  = eval_re_basis(re_coeff, Y_rot, erotc, ll)

b_size = size(eB_basis)[2]

difference_vec = [[Q*eB_basis[:,i] - eB_basis_rot[:,i]] for i in 1:b_size]
size(difference_vec)

print("Rotation symmetric: ", all([norm(difference_vec[i]) < tol for i in 1:b_size]),"\n")
n  = length(Y)

p = StatsBase.sample(1:4,4, replace=false)
eB_pi_basis  = eval_re_basis(re_coeff, Y, erotc, ll)
eB_pi_basis_rot  = eval_re_basis(re_coeff, Y_rot, erotc, ll)

difference_vec_perm = [[Q*eB_basis[:,i] - eB_basis_rot[:,i]] for i in 1:b_size]

[ (eval_re_basis(re_coeff, Y[p], erotc, ll) - eB_basis)  for p in Combinatorics.permutations(1:n)]

print("Rotation symmetric: ", all([norm(difference_vec_perm)< tol for i in 1:b_size]))

sample([rng], a, [wv::AbstractWeights], dims::Dims; replace=true, ordered=false)

end

difference_vec_perm = [Y[p]  for p in Combinatorics.permutations(1:n)]


text:
module ClickHouse
using Dates
using CategoricalArrays
include("Defines.jl")
include("Exceptions.jl")
include("tcp/tcp.jl")
include("columns/columns.jl")
include("Connect.jl")
include("Query.jl")

export ClickHouseSock
export Block
export select
export select_callback
export select_channel
export select_df
export insert
export execute
export connect, connect!
export ping
export ClickHouseServerException
export is_connected
export is_busy

end # module


text:
# io/z-list.jl

include("caller_name.jl")
include("fld-read.jl")
include("fld-write.jl")
include("ge-mri/pfile.jl")
#include("ge-mri/rdb-26_002.jl") # not
include("ge-mri/read_rdb_hdr.jl")
include("ir_dump.jl")
include("prompt.jl")
include("shows.jl")


text:
import Clustering, NearestNeighbors, LinearAlgebra
import Statistics
import StaticArrays

function get_transform_between_two_registered_clouds(c1, c2)
    centroid_1 = mean(c1, dims=2)
    centroid_2 = mean(c2, dims=2)
    centered_1 = c1 .- centroid_1
    centered_2 = c2 .- centroid_2
    H = centered_1 * centered_2'
    s = LinearAlgebra.svd(H)
    R = (s.V * permutedims(s.U))
    if LinearAlgebra.det(R) < 0
        s = LinearAlgebra.svd(R)
        V = s.V
        V[:,3] .= -1.0 .* V[:,3]
        R = V * permutedims(s.U)
    end
    rot = Rotations.RotMatrix(SMatrix{3,3}(R))

    T = (centroid_2 .- rot * centroid_1)
    T = T[:]
    Pose(T, rot)
end

"""
Returns the rigid transformation `T::Pose` such that `T * csource` is as close
as possible to `ctarget`.  Here `ctarget` and `csource` are point clouds, that
is, 3×something matrices where each column represents a point.
"""
function icp(c1, c2; iterations=10, c1_tree=nothing)
    T = IDENTITY_POSE
    if isnothing(c1_tree)
        c1_tree = NearestNeighbors.KDTree(c1)
    else
        @assert length(c1_tree.data) == size(c1)[2]
    end
    for _ in 1:iterations
        idxs, dists = NearestNeighbors.nn(c1_tree, c2)
        neighbors = c1[:,idxs]
        new_T = get_transform_between_two_registered_clouds(c2, neighbors)
        c2 = new_T.pos .+ (new_T.orientation * c2)
        T = get_c_relative_to_a(new_T, T)
    end

    T
end

function icp_object_pose(init_pose, obs, get_cloud_func;
    c1_tree=nothing, random_rotation_initialization=false, outer_iterations=5, iterations=10)
    p = init_pose
    if random_rotation_initialization
        p = Pose(p.pos, uniform_rot3())
    end

    if isnothing(c1_tree)
        c1_tree = NearestNeighbors.KDTree(obs)
    else
        @assert length(c1_tree.data) == size(obs)[2]
    end

    c = get_cloud_func(p)
    for _ in 1:outer_iterations
        T = icp(obs, c; iterations=iterations, c1_tree=c1_tree)
        p = get_c_relative_to_a(T,p)
        c = get_cloud_func(p)
    end

    p
end

function icp_project_to_planar_contact(
    obs, get_cloud_func, parent_pose, parent_box, parent_face, child_box, child_face,
    x,y,ang; c1_tree=nothing, outer_iterations=5, iterations=10)

    if isnothing(c1_tree)
        c1_tree = NearestNeighbors.KDTree(obs)
    else
        @assert length(c1_tree.data) == size(obs)[2]
    end

    contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
    p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
    c = get_cloud_func(p)

    for _ in 1:outer_iterations
        T = icp(obs, c; iterations=iterations, c1_tree=c1_tree)
        p = get_c_relative_to_a(T,p)

        (parent_face_idx, child_face_idx),_ = get_closest_planar_contact(parent_pose, parent_box, p, child_box)
        new_parent_face, new_child_face = S.BOX_SURFACE_IDS[parent_face_idx], S.BOX_SURFACE_IDS[child_face_idx]
        pc = S.closestApproximatingContact(
            (parent_pose * S.getContactPlane(parent_box, new_parent_face)) \
            (p * S.getContactPlane(child_box, new_child_face)),
        )
        x,y,ang = pc.x, pc.y, pc.angle
        contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang #= And no slack term. =#))
        p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)

        c = get_cloud_func(p)
    end

    (x,y,ang)
end

export get_transform_between_two_registered_clouds, icp, icp_object_pose

PyCall.py"""
import cv2 as cv
import numpy as np
def registerModelToScene(c1, c2):
    icp = cv.ppf_match_3d_ICP(6)
    return icp.registerModelToScene(np.array(c1).transpose().astype(np.float32),
                             np.array(c2).transpose().astype(np.float32))
def registerModelToSceneWithPoses(c1, c2, poses):
    icp = cv.ppf_match_3d_ICP(100)
    return icp.registerModelToScene(np.array(c1).transpose().astype(np.float32),
                             np.array(c2).transpose().astype(np.float32),
                             poses)
def pose_R_t_to_cv2_pose(R,t):
    p = cv.ppf_match_3d_Pose3D()
    p.updatePose(np.array(R), np.array(t))
    return p
"""
function cv2_pose_to_pose(pose::PyCall.PyObject)
   mat = pose.pose
   Pose(mat[:,4], RotMatrix3{Float64}(mat[1:3,1:3]))
end
function cv2_pose_to_pose(pose::Matrix)
   mat = pose
   Pose(mat[:,4], RotMatrix3{Float64}(mat[1:3,1:3]))
end
function pose_to_cv2_pose(pose)
    PyCall.py"pose_R_t_to_cv2_pose"(pose.orientation, pose.pos)
end

export icp_tree, icp_object_poses


text:
export rational_reconstruction, farey_lift, berlekamp_massey

## given some r/s = a mod b and deg(r) = n, deg(s) <= m find r,s
## a and b better be polynomials in the same poly ring.
## seems to work for Q (Qx) and Fp experimentally
#
# possibly should be rewritten to use an asymptotically (and practically)
# faster algorithm. For Q possibly using CRT and fast Fp techniques
# Algorithm copies from the bad-primes paper

@doc Markdown.doc"""
    rational_reconstruction(a::PolyElem{S}, b::PolyElem{S}, n::Int, m::Int)

 Returns true and x, y s.th. ay = x mod b and degree(x) <= n, degree(y) <= m
   or false (and garbage) if this is not possible.
"""
function rational_reconstruction(a::PolyElem{S}, b::PolyElem{S}, n::Int, m::Int) where S
  R = a.parent
  if degree(a) <= n return true, a, R(1); end

  M = zero_matrix(R, 2, 2)
  M[1,1] = b
  M[2,1] = a
  M[2,2] = R(1)

  T = zero_matrix(R, 2, 2)
  T[1,2] = R(1)
  T[2,1] = R(1)
  while degree(M[2,1]) > n
    q, r = divrem(M[1,1], M[2,1])
    T[2,2] = -q
    M = T*M
  end
  if degree(M[2,2]) <= m
    return true, M[2,1], M[2,2]
  end

  return false, M[2,1], M[2,2]
end

@doc Markdown.doc"""
  rational_reconstruction{S}(a::PolyElem{S}, b::PolyElem{S})

 Returns true and x/y s.th. ay = x mod b and degree(x), degree(y) <= degree(b)/2
   or false (and garbage) if this is not possible. Shortcut to the more general function.
"""
function rational_reconstruction(a::PolyElem{T}, b::PolyElem{T}) where T
  return rational_reconstruction_subres(a, b)
end
function rational_reconstruction(a::fmpq_poly, b::fmpq_poly)
  return rational_reconstruction_mod(a, b)
end

#Note: this is not the best (fastest) algorithm, not the most general
#      signature, not the best (fastest) interface, ....
#However: for now it works.

@doc Markdown.doc"""
    rational_reconstruction(a::fmpz, b::fmpz)
    rational_reconstruction(a::Integer, b::Integer)

Tries to solve ay=x mod b for x,y < sqrt(M/2). If possible, returns
  (true, x, y) or (false, garbage) if not possible.
"""
function rational_reconstruction(a::fmpz, b::fmpz)
  res = fmpq()
  a = mod(a, b)
  fl = ccall((:fmpq_reconstruct_fmpz, libflint), Int, (Ref{fmpq}, Ref{fmpz}, Ref{fmpz}), res, a, b)
  return fl!=0, numerator(res), denominator(res)
end

function rational_reconstruction(a::Integer, b::Integer)
  return rational_reconstruction(fmpz(a), fmpz(b))
end

@doc Markdown.doc"""
    rational_reconstruction(a::fmpz, b::fmpz, N::fmpz, D::fmpz) -> Bool, fmpz, fmpz
Given $a$ modulo $b$ and $N>0$, $D>0$ such that $2ND<b$, find $|x|\le N$, $0<y\le D$
satisfying $x/y \equiv a \bmod b$ or $a \equiv ya \bmod b$.
"""
function rational_reconstruction(a::fmpz, b::fmpz, N::fmpz, D::fmpz)
  res = fmpq()
  a = mod(a, b)
  fl = ccall((:fmpq_reconstruct_fmpz_2, libflint), Int, (Ref{fmpq}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), res, a, b, N, D)
  return fl!=0, numerator(res), denominator(res)
end

#Note: the vector version might be useful - or the mult be previous den version
#Note: missing reconstruction modulo a true ideal. W/o denominators

@doc Markdown.doc"""
    rational_reconstruction(a::nf_elem, b::fmpz)

Applies the rational_reconstruction function to each coefficient.
"""
function rational_reconstruction(a::nf_elem, b::fmpz)
  K= parent(a)
  Qx = parent(K.pol)
  res = Qx(0)
  for i=0:degree(K)-1
    c = coeff(a, i)
    if iszero(c)
      continue
    end
    @assert denominator(c) == 1
    fl, x, y = rational_reconstruction(numerator(c), b)
    if !fl
      return false, K(res)
    end
    setcoeff!(res, i, x//y)
  end
  return true, K(res)
end

# to appease the Singular crowd...
farey_lift = rational_reconstruction

# in at least 2 examples produces the same result as Magma
# can do experiments to see if dedicated Berlekamp Massey would be
# faster as well as experiments if Berlekamp Massey yields faster
# rational_reconstruction as well.
# Idea of using the same agorithm due to E. Thome
#

function berlekamp_massey_recon(a::Array{T, 1}) where T
  Rx,x = PolynomialRing(parent(a[1]), cached=false)
  f = Rx(a)
  xn= x^length(a)

  fl, n, d = rational_reconstruction(f, xn)
  if fl
    return true, d*(inv(trailing_coefficient(d)))
  else
    return false, Rx(0)
  end
end

###############################################################################
#                 univariate farey lift Michael Monagan algorithm             #
# from Dereje
###############################################################################

function rational_reconstruction_subres(g::PolyElem{T}, f::PolyElem{T}, bnd::Int = -1) where T
    # the denominator is normalized
    R_2 = g.parent
    r_1 = R_2(1); t_1 = R_2(0)
    r_m = r_1;t_m = r_1
    q1 = t_1; q_m = r_1

    if g==t_1
       return true, t_1,r_1
    end

    if 2*degree(g) < degree(f)
        return true, g, r_1
    end

    N1 = R_2(inv(lead(g))); r2 = g*N1
    r1 = f* R_2(inv(lead(f))); t1 = t_1;
    t2 = N1; i = 0
    l_rt = []
    deg_f = degree(f)
    while r2!=0
        i += 1
        q1,r = divrem(r1,r2); r3=r2
        if r==0
           N1 = R_2(1)
        else
           N1 = R_2(inv(lead(r)))
        end
        r2=r*N1; r1=r3
        r3=t2; t2=(t1-q1*t2)*N1; t1=r3

        if bnd>-1 && degree(r1) == bnd
          return true, r1, t1
        end

        if i == div(deg_f,2)
           l_rt = [r1, t1]
        end
        if degree(q1) > degree(q_m)
           (q_m, r_m, t_m) = (q1,r1,t1)
        end
    end

    if(degree(q_m)==1)
         if gcd(l_rt[1], l_rt[2])!=1
              return false, l_rt[1], l_rt[2]
         else
              return  true, l_rt[1], l_rt[2]
         end
    else
        if gcd(r_m, t_m) == 1
           return true, r_m, t_m
        else
           return false, r_m, t_m
        end
    end
end
###############################################################################
#                 modular univariate farey lift                               #
###############################################################################

function rational_reconstruction_mod(g::fmpq_poly, f::fmpq_poly)
  p = next_prime(fmpz(p_start))
  n, p = _inner_modp_results(g, f, p)  # mainly used to find the correct
                                       # bound n and a starting p
  kp = 10
  L =[]
  pp = FlintZZ(1)
  j = 0
  while true
    kp = 2*kp
    L = _modp_results(g,f,p,kp, n)
    p = L[4]
    if j==0
       N = L[1]; D = L[2]; pp = L[3]
    else
       N,_ = induce_crt(N, pp, L[1], L[3])
       D,pp = induce_crt(D, pp, L[2], L[3])
    end
    j=1
    fl, nu_rat_f = induce_rational_reconstruction(N, FlintZZ(pp))
    if fl
      fl, de_rat_f = induce_rational_reconstruction(D, FlintZZ(pp))
      if fl
        if (de_rat_f*g) % f == nu_rat_f
          return true,  nu_rat_f, de_rat_f
        end
      end
    end
    p = next_prime(p)
  end
end

################################################################################

function _modp_results(g::fmpq_poly,f::fmpq_poly, p::fmpz, M::Int, n::Int)
   l1 = nmod_poly[]; l2 = nmod_poly[];l3 = fmpz[]
   L = listprimes([f,g], p, M)
   for j in 1:length(L)
     Rp, t = PolynomialRing(ResidueRing(FlintZZ, L[j], cached=false), cached=false)
     gp = Rp(g)
     fp = Rp(f)
     fl, nu_p, de_p = rational_reconstruction_subres(gp, fp, n)
     if fl
        ut = Rp(inv(lead(de_p)))
        push!(l1, ut*nu_p)
        push!(l2, ut*de_p)
        push!(l3, L[j])
     end
   end
   c = crt_env(l3)
   nu = induce_crt(l1, c)
   du = induce_crt(l2, c)
   return nu, du, c.pr[end], L[end]
end

function _inner_modp_results(g::fmpq_poly,f::fmpq_poly, p::fmpz)
   while true
     if testPrime_jl(f,p) == true && testPrime_jl(g,p) == true
         Rp, t = PolynomialRing(ResidueRing(FlintZZ, p, cached=false), cached=false)
         gp = Rp(g)
         fp = Rp(f)
         fl, nu_p, de_p = rational_reconstruction_subres(gp, fp)
         if fl
             return degree(nu_p), p
         end
     end
     p = next_prime(p)
   end
end

###############################################################################

function berlekamp_massey(L::Array{T, 1}) where T
  return berlekamp_massey_naive(L)
end
function berlekamp_massey(L::Array{fmpq, 1})
  return berlekamp_massey_mod(L)
end

################################################################################
#                         Berlekamp Massey Algorithm                           #
################################################################################
function berlekamp_massey_naive(L::Array{T, 1}) where T
     R_s = parent(L[1])
     lg = length(L)
     L = [R_s(L[lg-i]) for i in 0:lg-1]
     Ry, Y = PolynomialRing(R_s, "Y", cached=false)
     g = Ry(L)
     if iszero(g)
       return true, g
     end
     f = Y^lg
     N = R_s(inv(lead(g))); g1 = g*N
     v0 = Ry(); v1 = Ry(1)
     while lg <= 2*degree(g1)
       q,r = divrem(f,g1)
       if r==0
          N = R_s(1)
       else
          N = R_s(inv(lead(r)))
       end
       v = (v0-q*v1)*N
       v0 = v1; v1 = v; f = g1; g1= r*N
     end
     return true, divexact(v1, lead(v1))
end

###############################################################################
#                 modular Berlekamp algorithm                                 #
###############################################################################

function berlekamp_massey_mod(L::Array{fmpq, 1})
  Rf = parent(L[1])
#  L = [Rf(L[i]) for i in 1:length(L)]
  Rc, Y = PolynomialRing(Rf, "Y", cached=false)
  f = Rc(L)
  if iszero(f)
    return true, f
  end
  p = next_prime(fmpz(p_start))
  kp = 10
  pp = FlintZZ(1)
  j = 0
  local N
  while true
    kp = 2*kp
    L = _modpResults(f,p,kp)
    p = L[3]
    if j==0
       N = L[1]; pp = L[2]
    else
       N, pp = induce_crt(N, pp, L[1], L[2])
      j=1
    end
    fl, nu_rat_f = induce_rational_reconstruction(N, FlintZZ(pp))
    if fl
      return true, nu_rat_f
      #the check for roots is ONLY useful in multivariate interpolation
      #in general, the poly can be anyhting of course
      if length(factor(nu_rat_f)) == degree(nu_rat_f) #TODO write and use roots
          return true,  nu_rat_f
      end
    end
    p = next_prime(p)
  end
end

################################################################################

function _modpResults(f, p::fmpz, M::Int)

   Rc = f.parent
   l1 = fmpz_mod_poly[]; l3 = fmpz[]
   Np = listprimes([f], p, M)
   Zx, Y = PolynomialRing(FlintZZ, "Y", cached=false)
   for j in 1:length(Np)
     RNp = ResidueRing(FlintZZ, Np[j], cached=false)
     Rp, t = PolynomialRing(RNp, "t", cached=false)
     fp = Rp(f)
     L1 = Nemo.fmpz_mod[]
     for i in 0:degree(fp)
        push!(L1, coeff(fp, i))
     end
     ff = berlekamp_massey_naive(L1)
     push!(l1, ff)
     push!(l3, Np[j])
   end
   c = crt_env(l3)
   nu = induce_crt(l1, c)
   return nu, c.pr[end], Np[end]
end

################################################################################

function testPrime_jl(f::fmpq_poly, p::fmpz)
    # BAD!!! missing: num_coeff(f, i)
    nd = denominator(f)
    fg = nd*f
    return !(divides(nd, p)[1]) || !(divides(numerator(lead(fg)), p)[1])
end

################################################################################

function listprimes(f::Array{fmpq_poly, 1}, p::fmpz, M::Int)
   # static
   i=0; L = fmpz[]
   while true
    i += 1
    if all(x-> testPrime_jl(x,p), f)
       push!(L, p)
    end
    if i == M
      return L
    end
    p = next_prime(p)
  end
end

################################################################################

function induce_crt(L::Array{nmod_poly, 1}, c::crt_env{fmpz})
  Zx, x = FlintZZ["x"]
  res = Zx()
  m = maximum(degree(x) for x = L)

  for i=0:m
    setcoeff!(res, i, crt([lift(coeff(x, i)) for x =L], c))
  end
  return res
end