text:
using Test
# Write a test that compares the momintegrals_nested approach with applying the
# Wilton singularity extraction to a non-conforming mesh.

# Conclusion: WiltonSEStrategy should not be applied to deal with overlapping
# charts in the case of geometrically non-conforming meshes; the result of the
# inner integral is not analytic, crippling the accuracy of the outer integration,
# which is done with triangular Gauss-Legendre points.

using CompScienceMeshes, BEAST
using LinearAlgebra

o = point(0,0,0)
c = point(1/3,1/3,0)
e = point(1/2,0,0)
d = point(1/3,-1/3,0)

# Relative weights chosen to make the two terms equally important
𝒜 = Maxwell3D.singlelayer(wavenumber=0.0, alpha=22.0, beta=1.0)
𝒳 = BEAST.RTRefSpace{Float64}()
T = scalartype(𝒜, 𝒳, 𝒳)

test_chart = simplex(e,x̂,c)
trial_chart = simplex(o,x̂,ŷ)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

function BEAST.quaddata(op::BEAST.MaxwellOperator3D,
    test_local_space::BEAST.RefSpace, trial_local_space::BEAST.RefSpace,
    test_charts, trial_charts)

    a, b = 0.0, 1.0
    tqd = quadpoints(test_local_space, test_charts, (12,12))
    bqd = quadpoints(trial_local_space, trial_charts, (13,13))
    leg = (BEAST._legendre(10,a,b), BEAST._legendre(10,a,b), BEAST._legendre(10,a,b),)
    return (tpoints=tqd, bpoints=bqd, gausslegendre=leg)
end

sauterschwab = BEAST.SauterSchwabQuadrature.CommonFace(nothing)
out_ss = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals_nested!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals_nested!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw,wiltonsingext)

@show norm(out_ss-out_dw) / norm(out_dw)

test_chart = simplex(e,x̂,c)
trial_chart = simplex(e,x̂,c)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonFace(BEAST._legendre(10,0.0,1.0))
out_ss1 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss1,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw1 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw1,wiltonsingext)

@test out_ss1 ≈ out_dw1 rtol=3e-6

# Test accuracy CommonEdge
test_chart = simplex(e,x̂,c)
trial_chart = simplex(c,x̂,(x̂+ŷ)/2)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0))
out_ss2 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss2,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw2 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw2,wiltonsingext)

@test out_ss2 ≈ out_dw2 rtol=4e-6

test_chart = simplex(e,x̂,c)
trial_chart = simplex(o,e,c)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0))
out_ss3 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss3,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw3 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw3,wiltonsingext)

@test out_ss3 ≈ out_dw3 rtol=3e-6

# Test the CommonVertex cases

test_chart = simplex(e,x̂,c)
trial_chart = simplex(c,(x̂+ŷ)/2,ŷ)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
out_ss4 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss4,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw4 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw4,wiltonsingext)

@test out_ss4 ≈ out_dw4 rtol=2e-9

test_chart = simplex(e,x̂,c)
trial_chart = simplex(c,ŷ,ŷ/2)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
out_ss5 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss5,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw5 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw5,wiltonsingext)

@test out_ss4 ≈ out_dw4 rtol=3e-8

test_chart = simplex(e,x̂,c)
trial_chart = simplex(c,ŷ/2,o)

test_quadpoints  = BEAST.quadpoints(𝒳, [test_chart],  (12,))[1,1]
trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]

sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
out_ss6 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss6,sauterschwab)

wiltonsingext = BEAST.WiltonSEStrategy(test_quadpoints, BEAST.DoubleQuadStrategy(test_quadpoints, trial_quadpoints))
out_dw6 = zeros(T, numfunctions(𝒳), numfunctions(𝒳))
BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw6,wiltonsingext)

@test out_ss4 ≈ out_dw4 rtol=6e-8

out_ss = out_ss1 + out_ss2 + out_ss3 + out_ss4 + out_ss5 + out_ss6
out_dw = out_dw1 + out_dw2 + out_dw3 + out_dw4 + out_dw5 + out_dw6

@test out_ss ≈ out_dw rtol=4e-7


text:
# Test whether degeneracy increases with numgates (numinteriors).
# Test whether degeneracy increases with sporadic mutation of the goal.
export recover_phenotype

function test_degen( numinputs::Int64, numoutputs::Int64, numints_rng::IntRange,
    numtrials::Int64, maxsteps::Int64 )
  df = DataFrame()
  df.numinputs = Int64[]
  df.numoutputs = Int64[]
  df.numints = Int64[]
  df.levsback = Int64[]
  df.degen = Float64[]
  for numints in numints_rng
    levsback = Int(ceil(numints/2))
    sum_degen = 0.0
    for t = 1:numtrials
      p = Parameters( numinputs, numoutputs, numints, levsback )
      c = random_chromosome(p)
      g = randgoal(numinputs,numoutputs)
      new_c = neutral_evolution( c, g, maxsteps )
      sum_degen += degeneracy( c )
    end
    row = (numinputs,numoutputs,numints,levsback,sum_degen/numtrials)
    push!(df,row)
  end
  df
end

# Helper function for test_degen()
function run_test_degen(numinputs::Int64, numoutputs::Int64, numints::Int64, maxsteps::Int64 )
  p = Parameters( numinputs, numoutputs, numints, levsback )
  c = random_chromosome(p)
  g = randgoal(numinputs,numoutputs)
  new_c = neutral_evolution( c, g, maxsteps )
  sum_degen += degeneracy( c )
end

# Starts with a chromsome c.
# For each trial, mutate the chromosome, and use neutral evolution to try evolve the mutated chromosome
#   to output the phenotype of c.  Up to maxtries attempts at neutral evolution are made
# Results the are the sum of the evolutionary steps made by neutral evolution
#   and the sum of the number of tries at neutral evolution.
function recover_phenotype( c::Chromosome, maxsteps::Int64, maxtrials::Int64, maxtries::Int64 )
  p = c.params
  funcs = default_funcs(p.numinputs)
  g = output_values(c)
  sumsteps = 0
  sumtries = 0
  for t = 1:maxtrials
    (mc,active) = mutate_chromosome!(deepcopy(c),funcs)
    (nc,steps) = neutral_evolution( mc, g, maxsteps )
    sumsteps += steps
    tries = 1
    while tries < maxtries && steps == maxsteps
      (nc,steps) = neutral_evolution( mc, g, maxsteps )
      sumsteps += steps
      tries += 1
    end
    if tries == maxtries
      println("tries: ",tries,"  maxtries: ",maxtries)
    end
    sumtries += tries
  end
  return (sumsteps,sumtries)
end

#=
# Moved to Evolve.jl.
function lambda_evolution( c::Chromosome, g::Goal, maxsteps::Integer, mutrate::Float64 )
  p = c.params
  p.mutrate = mutrate
  funcs = default_funcs(p.numinputs)
  poisson_lambda = p.mutrate*p.numinteriors
  X = Poisson( poisson_lambda )
  step = 0
  ov = output_values( c)
  current_distance = hamming_distance( ov, g, c.params.numinputs )
  #print("ov: ",ov,"  cur_dist: ",current_distance,"  "); print_circuit(c)
  while step < maxsteps && ov != g
    chrome_list = Chromosome[]
    dist_list = Float64[]
    for i = 1:p.lambda
      step += 1
      num_mutations =  rand(X)
      new_c = deepcopy(c)
      for m = 1:num_mutations
        mutate_chromosome!( new_c, funcs )
      end
      new_ov = output_values( new_c)
      new_dist = hamming_distance( new_ov, g, c.params.numinputs )
      #print("new_ov: ",new_ov,"  new_dist: ",new_dist,"  "); print_circuit(new_c)
      push!(chrome_list,new_c)
      push!(dist_list,new_dist)
    end
    (best_dist,ind) = findmin( dist_list )
    if best_dist <= current_distance
      c = chrome_list[ind]
      ov = output_values( c )
      current_distance = hamming_distance( ov, g, c.params.numinputs )
      #print("ov: ",ov,"  cur_dist: ",current_distance,"  "); print_circuit(c)
    end
  end # while
  if step == maxsteps
    println("lambda evolution failed with ",step," steps for goal: ",g)
    return (c, step)
  else
    println("lambda evolution succeeded at step ",step," for goal: ",g)
    @assert output_values(c) == g
    return (c, step)
  end
end
=#

function compare_lambda_neutral( p::Parameters, g::Goal, trials::Int64, maxsteps::Int64, mutrate_range::CGP.FloatRange )
  @assert length(collect(mutrate_range)) == 3
  df = DataFrame()
  df.numinputs = Int64[]
  df.numoutputs = Int64[]
  df.numinteriorss = Int64[]
  df.numlevelsback = Int64[]
  df.mutrate = Float64[]
  df.neutral_fails = Int64[]
  df.neutral_steps= Int64[]
  df.lamda1_fails = Int64[]
  df.lamda1_steps = Int64[]
  df.lamda2_fails = Int64[]
  df.lamda2_steps = Int64[]
  df.lamda3_fails = Int64[]
  df.lamda3_steps = Int64[]
  for t = 1:numtrials
    c = random_chromosome(p)
    (neut_c,neut_steps) = neutral_evolution( deepcopy(c), g, maxsteps )
    mr = collect(mutrate_range)
    (lambda_c1,lambda_steps1) = lambda_evolution( deepcopy(c), g, maxsteps, mr[1] )
    (lambda_c2,lambda_steps2) = lambda_evolution( deepcopy(c), g, maxsteps, mr[2] )
    (lambda_c3,lambda_steps3) = lambda_evolution( deepcopy(c), g, maxsteps, mr[3] )
    df_row = (
      p.numinputs,
      p.numoutputs
    )
  end
end

ext: jl


text:
module PLaplacianWithAutodiffTests

using Test
using Gridap

u(x) = x[1] + x[2]
f(x) = 0
const p = 3
@law j(∇u) = norm(∇u)^(p-2) * ∇u

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 3
T = Float64

V = TestFESpace(
 model=model,
 order=order,
 reffe=:Lagrangian,
 conformity=:H1,
 valuetype=T,
 dirichlet_tags="boundary")

U = TrialFESpace(V,u)

degree = 2*order
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

res(u,v) = ∇(v)⋅j(∇(u)) - v*f

t_Ω = FETerm(res,trian,quad)

op = FEOperator(U,V,t_Ω)

nls = NLSolver(show_trace=false, method=:newton)
solver = FESolver(nls)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)

e = u - uh

el2 = sqrt(sum(integrate(e*e,trian,quad)))
eh1 = sqrt(sum(integrate(e*e+∇(e)⋅∇(e),trian,quad)))

@test el2 < 1.e-8
@test eh1 < 1.e-7

end # module


text:
# This file is a part of NeighborhoodApproximationIndex.jl

using NeighborhoodApproximationIndex
using Test

using KCenters, SimilaritySearch, NeighborhoodApproximationIndex

#include("invindex.jl")
#include("testknr.jl")

@testset "NeighborhoodApproximationIndex.jl" begin
    dim = 4
    n = 100_000
    X = [randn(dim) for i in 1:n]
    Q = [randn(dim) for i in 1:100]
    dist = SqL2Distance()
    seq = ExhaustiveSearch(dist, X)
    kbuild, ksearch, t, recall = 3, 1, 1, 0.9
    P = Performance(seq, Q, 7)
    index = DeloneInvertedFile(dist, X; centers=:exhaustive, refs=:rand, k=kbuild, t=t)
    p = probe(P, copy(index, k=ksearch))
    @info "before optimization: $(index)" (recall=p.macrorecall, queries_per_second= 1 / p.searchtime, eval_ratio=p.evaluations / length(X))
    @test p.macrorecall > recall
    #optimize!(P, index, 0.9, verbose=true)
    #p = probe(P, index)

    #@info "after optimization: $(string(index))" (recall=p.macrorecall, queries_per_second= 1 / p.searchtime, eval_ratio=p.evaluations / length(X))
    # we can expect a small recall reduction since we are using a db's subset for optimizin the index

end


text:
function massmat20!(Out, In)

    Out[1+0] = 8.86049999999999;
    Out[1+1] = 0.15958366*cos(In[2]);
    Out[1+2] = 0.151604477*cos(In[3]);
    Out[1+3] = 0.143625294*cos(In[4]);
    Out[1+4] = 0.135646111*cos(In[5]);
    Out[1+5] = 0.127666928*cos(In[6]);
    Out[1+6] = 0.119687745*cos(In[7]);
    Out[1+7] = 0.111708562*cos(In[8]);
    Out[1+8] = 0.103729379*cos(In[9]);
    Out[1+9] = 0.095750196*cos(In[10]);
    Out[1+10] = 0.087771013*cos(In[11]);
    Out[1+11] = 0.07979183*cos(In[12]);
    Out[1+12] = 0.071812647*cos(In[13]);
    Out[1+13] = 0.063833464*cos(In[14]);
    Out[1+14] = 0.055854281*cos(In[15]);
    Out[1+15] = 0.047875098*cos(In[16]);
    Out[1+16] = 0.039895915*cos(In[17]);
    Out[1+17] = 0.031916732*cos(In[18]);
    Out[1+18] = 0.023937549*cos(In[19]);
    Out[1+19] = 0.015958366*cos(In[20]);
    Out[1+20] = 0.007979183*cos(In[21]);
    Out[1+21] = 0.15958366*cos(In[2]);
    Out[1+22] = 0.009399477574;
    Out[1+23] = 0.0089295036953*sin(In[2])*sin(In[3]) + 0.0089295036953*cos(In[2])*cos(In[3]);
    Out[1+24] = 0.0084595298166*sin(In[2])*sin(In[4]) + 0.0084595298166*cos(In[2])*cos(In[4]);
    Out[1+25] = 0.0079895559379*sin(In[2])*sin(In[5]) + 0.0079895559379*cos(In[2])*cos(In[5]);
    Out[1+26] = 0.0075195820592*sin(In[2])*sin(In[6]) + 0.0075195820592*cos(In[2])*cos(In[6]);
    Out[1+27] = 0.0070496081805*sin(In[2])*sin(In[7]) + 0.0070496081805*cos(In[2])*cos(In[7]);
    Out[1+28] = 0.0065796343018*sin(In[2])*sin(In[8]) + 0.0065796343018*cos(In[2])*cos(In[8]);
    Out[1+29] = 0.0061096604231*sin(In[2])*sin(In[9]) + 0.0061096604231*cos(In[2])*cos(In[9]);
    Out[1+30] = 0.0056396865444*sin(In[10])*sin(In[2]) + 0.0056396865444*cos(In[10])*cos(In[2]);
    Out[1+31] = 0.0051697126657*sin(In[11])*sin(In[2]) + 0.0051697126657*cos(In[11])*cos(In[2]);
    Out[1+32] = 0.004699738787*sin(In[12])*sin(In[2]) + 0.004699738787*cos(In[12])*cos(In[2]);
    Out[1+33] = 0.0042297649083*sin(In[13])*sin(In[2]) + 0.0042297649083*cos(In[13])*cos(In[2]);
    Out[1+34] = 0.0037597910296*sin(In[14])*sin(In[2]) + 0.0037597910296*cos(In[14])*cos(In[2]);
    Out[1+35] = 0.0032898171509*sin(In[15])*sin(In[2]) + 0.0032898171509*cos(In[15])*cos(In[2]);
    Out[1+36] = 0.0028198432722*sin(In[16])*sin(In[2]) + 0.0028198432722*cos(In[16])*cos(In[2]);
    Out[1+37] = 0.0023498693935*sin(In[17])*sin(In[2]) + 0.0023498693935*cos(In[17])*cos(In[2]);
    Out[1+38] = 0.0018798955148*sin(In[18])*sin(In[2]) + 0.0018798955148*cos(In[18])*cos(In[2]);
    Out[1+39] = 0.0014099216361*sin(In[19])*sin(In[2]) + 0.0014099216361*cos(In[19])*cos(In[2]);
    Out[1+40] = 0.0009399477574*sin(In[20])*sin(In[2]) + 0.0009399477574*cos(In[20])*cos(In[2]);
    Out[1+41] = 0.0004699738787*sin(In[21])*sin(In[2]) + 0.0004699738787*cos(In[21])*cos(In[2]);
    Out[1+42] = 0.151604477*cos(In[3]);
    Out[1+43] = 0.0089295036953*sin(In[2])*sin(In[3]) + 0.0089295036953*cos(In[2])*cos(In[3]);
    Out[1+44] = 0.0089295036953;
    Out[1+45] = 0.0084595298166*sin(In[3])*sin(In[4]) + 0.0084595298166*cos(In[3])*cos(In[4]);
    Out[1+46] = 0.0079895559379*sin(In[3])*sin(In[5]) + 0.0079895559379*cos(In[3])*cos(In[5]);
    Out[1+47] = 0.0075195820592*sin(In[3])*sin(In[6]) + 0.0075195820592*cos(In[3])*cos(In[6]);
    Out[1+48] = 0.0070496081805*sin(In[3])*sin(In[7]) + 0.0070496081805*cos(In[3])*cos(In[7]);
    Out[1+49] = 0.0065796343018*sin(In[3])*sin(In[8]) + 0.0065796343018*cos(In[3])*cos(In[8]);
    Out[1+50] = 0.0061096604231*sin(In[3])*sin(In[9]) + 0.0061096604231*cos(In[3])*cos(In[9]);
    Out[1+51] = 0.0056396865444*sin(In[10])*sin(In[3]) + 0.0056396865444*cos(In[10])*cos(In[3]);
    Out[1+52] = 0.0051697126657*sin(In[11])*sin(In[3]) + 0.0051697126657*cos(In[11])*cos(In[3]);
    Out[1+53] = 0.004699738787*sin(In[12])*sin(In[3]) + 0.004699738787*cos(In[12])*cos(In[3]);
    Out[1+54] = 0.0042297649083*sin(In[13])*sin(In[3]) + 0.0042297649083*cos(In[13])*cos(In[3]);
    Out[1+55] = 0.0037597910296*sin(In[14])*sin(In[3]) + 0.0037597910296*cos(In[14])*cos(In[3]);
    Out[1+56] = 0.0032898171509*sin(In[15])*sin(In[3]) + 0.0032898171509*cos(In[15])*cos(In[3]);
    Out[1+57] = 0.0028198432722*sin(In[16])*sin(In[3]) + 0.0028198432722*cos(In[16])*cos(In[3]);
    Out[1+58] = 0.0023498693935*sin(In[17])*sin(In[3]) + 0.0023498693935*cos(In[17])*cos(In[3]);
    Out[1+59] = 0.0018798955148*sin(In[18])*sin(In[3]) + 0.0018798955148*cos(In[18])*cos(In[3]);
    Out[1+60] = 0.0014099216361*sin(In[19])*sin(In[3]) + 0.0014099216361*cos(In[19])*cos(In[3]);
    Out[1+61] = 0.0009399477574*sin(In[20])*sin(In[3]) + 0.0009399477574*cos(In[20])*cos(In[3]);
    Out[1+62] = 0.0004699738787*sin(In[21])*sin(In[3]) + 0.0004699738787*cos(In[21])*cos(In[3]);
    Out[1+63] = 0.143625294*cos(In[4]);
    Out[1+64] = 0.0084595298166*sin(In[2])*sin(In[4]) + 0.0084595298166*cos(In[2])*cos(In[4]);
    Out[1+65] = 0.0084595298166*sin(In[3])*sin(In[4]) + 0.0084595298166*cos(In[3])*cos(In[4]);
    Out[1+66] = 0.0084595298166;
    Out[1+67] = 0.0079895559379*sin(In[4])*sin(In[5]) + 0.0079895559379*cos(In[4])*cos(In[5]);
    Out[1+68] = 0.0075195820592*sin(In[4])*sin(In[6]) + 0.0075195820592*cos(In[4])*cos(In[6]);
    Out[1+69] = 0.0070496081805*sin(In[4])*sin(In[7]) + 0.0070496081805*cos(In[4])*cos(In[7]);
    Out[1+70] = 0.0065796343018*sin(In[4])*sin(In[8]) + 0.0065796343018*cos(In[4])*cos(In[8]);
    Out[1+71] = 0.0061096604231*sin(In[4])*sin(In[9]) + 0.0061096604231*cos(In[4])*cos(In[9]);
    Out[1+72] = 0.0056396865444*sin(In[10])*sin(In[4]) + 0.0056396865444*cos(In[10])*cos(In[4]);
    Out[1+73] = 0.0051697126657*sin(In[11])*sin(In[4]) + 0.0051697126657*cos(In[11])*cos(In[4]);
    Out[1+74] = 0.004699738787*sin(In[12])*sin(In[4]) + 0.004699738787*cos(In[12])*cos(In[4]);
    Out[1+75] = 0.0042297649083*sin(In[13])*sin(In[4]) + 0.0042297649083*cos(In[13])*cos(In[4]);
    Out[1+76] = 0.0037597910296*sin(In[14])*sin(In[4]) + 0.0037597910296*cos(In[14])*cos(In[4]);
    Out[1+77] = 0.0032898171509*sin(In[15])*sin(In[4]) + 0.0032898171509*cos(In[15])*cos(In[4]);
    Out[1+78] = 0.0028198432722*sin(In[16])*sin(In[4]) + 0.0028198432722*cos(In[16])*cos(In[4]);
    Out[1+79] = 0.0023498693935*sin(In[17])*sin(In[4]) + 0.0023498693935*cos(In[17])*cos(In[4]);
    Out[1+80] = 0.0018798955148*sin(In[18])*sin(In[4]) + 0.0018798955148*cos(In[18])*cos(In[4]);
    Out[1+81] = 0.0014099216361*sin(In[19])*sin(In[4]) + 0.0014099216361*cos(In[19])*cos(In[4]);
    Out[1+82] = 0.0009399477574*sin(In[20])*sin(In[4]) + 0.0009399477574*cos(In[20])*cos(In[4]);
    Out[1+83] = 0.0004699738787*sin(In[21])*sin(In[4]) + 0.0004699738787*cos(In[21])*cos(In[4]);
    Out[1+84] = 0.135646111*cos(In[5]);
    Out[1+85] = 0.0079895559379*sin(In[2])*sin(In[5]) + 0.0079895559379*cos(In[2])*cos(In[5]);
    Out[1+86] = 0.0079895559379*sin(In[3])*sin(In[5]) + 0.0079895559379*cos(In[3])*cos(In[5]);
    Out[1+87] = 0.0079895559379*sin(In[4])*sin(In[5]) + 0.0079895559379*cos(In[4])*cos(In[5]);
    Out[1+88] = 0.0079895559379;
    Out[1+89] = 0.0075195820592*sin(In[5])*sin(In[6]) + 0.0075195820592*cos(In[5])*cos(In[6]);
    Out[1+90] = 0.0070496081805*sin(In[5])*sin(In[7]) + 0.0070496081805*cos(In[5])*cos(In[7]);
    Out[1+91] = 0.0065796343018*sin(In[5])*sin(In[8]) + 0.0065796343018*cos(In[5])*cos(In[8]);
    Out[1+92] = 0.0061096604231*sin(In[5])*sin(In[9]) + 0.0061096604231*cos(In[5])*cos(In[9]);
    Out[1+93] = 0.0056396865444*sin(In[10])*sin(In[5]) + 0.0056396865444*cos(In[10])*cos(In[5]);
    Out[1+94] = 0.0051697126657*sin(In[11])*sin(In[5]) + 0.0051697126657*cos(In[11])*cos(In[5]);
    Out[1+95] = 0.004699738787*sin(In[12])*sin(In[5]) + 0.004699738787*cos(In[12])*cos(In[5]);
    Out[1+96] = 0.0042297649083*sin(In[13])*sin(In[5]) + 0.0042297649083*cos(In[13])*cos(In[5]);
    Out[1+97] = 0.0037597910296*sin(In[14])*sin(In[5]) + 0.0037597910296*cos(In[14])*cos(In[5]);
    Out[1+98] = 0.0032898171509*sin(In[15])*sin(In[5]) + 0.0032898171509*cos(In[15])*cos(In[5]);
    Out[1+99] = 0.0028198432722*sin(In[16])*sin(In[5]) + 0.0028198432722*cos(In[16])*cos(In[5]);
    Out[1+100] = 0.0023498693935*sin(In[17])*sin(In[5]) + 0.0023498693935*cos(In[17])*cos(In[5]);
    Out[1+101] = 0.0018798955148*sin(In[18])*sin(In[5]) + 0.0018798955148*cos(In[18])*cos(In[5]);
    Out[1+102] = 0.0014099216361*sin(In[19])*sin(In[5]) + 0.0014099216361*cos(In[19])*cos(In[5]);
    Out[1+103] = 0.0009399477574*sin(In[20])*sin(In[5]) + 0.0009399477574*cos(In[20])*cos(In[5]);
    Out[1+104] = 0.0004699738787*sin(In[21])*sin(In[5]) + 0.0004699738787*cos(In[21])*cos(In[5]);
    Out[1+105] = 0.127666928*cos(In[6]);
    Out[1+106] = 0.0075195820592*sin(In[2])*sin(In[6]) + 0.0075195820592*cos(In[2])*cos(In[6]);
    Out[1+107] = 0.0075195820592*sin(In[3])*sin(In[6]) + 0.0075195820592*cos(In[3])*cos(In[6]);
    Out[1+108] = 0.0075195820592*sin(In[4])*sin(In[6]) + 0.0075195820592*cos(In[4])*cos(In[6]);
    Out[1+109] = 0.0075195820592*sin(In[5])*sin(In[6]) + 0.0075195820592*cos(In[5])*cos(In[6]);
    Out[1+110] = 0.0075195820592;
    Out[1+111] = 0.0070496081805*sin(In[6])*sin(In[7]) + 0.0070496081805*cos(In[6])*cos(In[7]);
    Out[1+112] = 0.0065796343018*sin(In[6])*sin(In[8]) + 0.0065796343018*cos(In[6])*cos(In[8]);
    Out[1+113] = 0.0061096604231*sin(In[6])*sin(In[9]) + 0.0061096604231*cos(In[6])*cos(In[9]);
    Out[1+114] = 0.0056396865444*sin(In[10])*sin(In[6]) + 0.0056396865444*cos(In[10])*cos(In[6]);
    Out[1+115] = 0.0051697126657*sin(In[11])*sin(In[6]) + 0.0051697126657*cos(In[11])*cos(In[6]);
    Out[1+116] = 0.004699738787*sin(In[12])*sin(In[6]) + 0.004699738787*cos(In[12])*cos(In[6]);
    Out[1+117] = 0.0042297649083*sin(In[13])*sin(In[6]) + 0.0042297649083*cos(In[13])*cos(In[6]);
    Out[1+118] = 0.0037597910296*sin(In[14])*sin(In[6]) + 0.0037597910296*cos(In[14])*cos(In[6]);
    Out[1+119] = 0.0032898171509*sin(In[15])*sin(In[6]) + 0.0032898171509*cos(In[15])*cos(In[6]);
    Out[1+120] = 0.0028198432722*sin(In[16])*sin(In[6]) + 0.0028198432722*cos(In[16])*cos(In[6]);
    Out[1+121] = 0.0023498693935*sin(In[17])*sin(In[6]) + 0.0023498693935*cos(In[17])*cos(In[6]);
    Out[1+122] = 0.0018798955148*sin(In[18])*sin(In[6]) + 0.0018798955148*cos(In[18])*cos(In[6]);
    Out[1+123] = 0.0014099216361*sin(In[19])*sin(In[6]) + 0.0014099216361*cos(In[19])*cos(In[6]);
    Out[1+124] = 0.0009399477574*sin(In[20])*sin(In[6]) + 0.0009399477574*cos(In[20])*cos(In[6]);
    Out[1+125] = 0.0004699738787*sin(In[21])*sin(In[6]) + 0.0004699738787*cos(In[21])*cos(In[6]);
    Out[1+126] = 0.119687745*cos(In[7]);
    Out[1+127] = 0.0070496081805*sin(In[2])*sin(In[7]) + 0.0070496081805*cos(In[2])*cos(In[7]);
    Out[1+128] = 0.0070496081805*sin(In[3])*sin(In[7]) + 0.0070496081805*cos(In[3])*cos(In[7]);
    Out[1+129] = 0.0070496081805*sin(In[4])*sin(In[7]) + 0.0070496081805*cos(In[4])*cos(In[7]);
    Out[1+130] = 0.0070496081805*sin(In[5])*sin(In[7]) + 0.0070496081805*cos(In[5])*cos(In[7]);
    Out[1+131] = 0.0070496081805*sin(In[6])*sin(In[7]) + 0.0070496081805*cos(In[6])*cos(In[7]);
    Out[1+132] = 0.0070496081805;
    Out[1+133] = 0.0065796343018*sin(In[7])*sin(In[8]) + 0.0065796343018*cos(In[7])*cos(In[8]);
    Out[1+134] = 0.0061096604231*sin(In[7])*sin(In[9]) + 0.0061096604231*cos(In[7])*cos(In[9]);
    Out[1+135] = 0.0056396865444*sin(In[10])*sin(In[7]) + 0.0056396865444*cos(In[10])*cos(In[7]);
    Out[1+136] = 0.0051697126657*sin(In[11])*sin(In[7]) + 0.0051697126657*cos(In[11])*cos(In[7]);
    Out[1+137] = 0.004699738787*sin(In[12])*sin(In[7]) + 0.004699738787*cos(In[12])*cos(In[7]);
    Out[1+138] = 0.0042297649083*sin(In[13])*sin(In[7]) + 0.0042297649083*cos(In[13])*cos(In[7]);
    Out[1+139] = 0.0037597910296*sin(In[14])*sin(In[7]) + 0.0037597910296*cos(In[14])*cos(In[7]);
    Out[1+140] = 0.0032898171509*sin(In[15])*sin(In[7]) + 0.0032898171509*cos(In[15])*cos(In[7]);
    Out[1+141] = 0.0028198432722*sin(In[16])*sin(In[7]) + 0.0028198432722*cos(In[16])*cos(In[7]);
    Out[1+142] = 0.0023498693935*sin(In[17])*sin(In[7]) + 0.0023498693935*cos(In[17])*cos(In[7]);
    Out[1+143] = 0.0018798955148*sin(In[18])*sin(In[7]) + 0.0018798955148*cos(In[18])*cos(In[7]);
    Out[1+144] = 0.0014099216361*sin(In[19])*sin(In[7]) + 0.0014099216361*cos(In[19])*cos(In[7]);
    Out[1+145] = 0.0009399477574*sin(In[20])*sin(In[7]) + 0.0009399477574*cos(In[20])*cos(In[7]);
    Out[1+146] = 0.0004699738787*sin(In[21])*sin(In[7]) + 0.0004699738787*cos(In[21])*cos(In[7]);
    Out[1+147] = 0.111708562*cos(In[8]);
    Out[1+148] = 0.0065796343018*sin(In[2])*sin(In[8]) + 0.0065796343018*cos(In[2])*cos(In[8]);
    Out[1+149] = 0.0065796343018*sin(In[3])*sin(In[8]) + 0.0065796343018*cos(In[3])*cos(In[8]);
    Out[1+150] = 0.0065796343018*sin(In[4])*sin(In[8]) + 0.0065796343018*cos(In[4])*cos(In[8]);
    Out[1+151] = 0.0065796343018*sin(In[5])*sin(In[8]) + 0.0065796343018*cos(In[5])*cos(In[8]);
    Out[1+152] = 0.0065796343018*sin(In[6])*sin(In[8]) + 0.0065796343018*cos(In[6])*cos(In[8]);
    Out[1+153] = 0.0065796343018*sin(In[7])*sin(In[8]) + 0.0065796343018*cos(In[7])*cos(In[8]);
    Out[1+154] = 0.0065796343018;
    Out[1+155] = 0.0061096604231*sin(In[8])*sin(In[9]) + 0.0061096604231*cos(In[8])*cos(In[9]);
    Out[1+156] = 0.0056396865444*sin(In[10])*sin(In[8]) + 0.0056396865444*cos(In[10])*cos(In[8]);
    Out[1+157] = 0.0051697126657*sin(In[11])*sin(In[8]) + 0.0051697126657*cos(In[11])*cos(In[8]);
    Out[1+158] = 0.004699738787*sin(In[12])*sin(In[8]) + 0.004699738787*cos(In[12])*cos(In[8]);
    Out[1+159] = 0.0042297649083*sin(In[13])*sin(In[8]) + 0.0042297649083*cos(In[13])*cos(In[8]);
    Out[1+160] = 0.0037597910296*sin(In[14])*sin(In[8]) + 0.0037597910296*cos(In[14])*cos(In[8]);
    Out[1+161] = 0.0032898171509*sin(In[15])*sin(In[8]) + 0.0032898171509*cos(In[15])*cos(In[8]);
    Out[1+162] = 0.0028198432722*sin(In[16])*sin(In[8]) + 0.0028198432722*cos(In[16])*cos(In[8]);
    Out[1+163] = 0.0023498693935*sin(In[17])*sin(In[8]) + 0.0023498693935*cos(In[17])*cos(In[8]);
    Out[1+164] = 0.0018798955148*sin(In[18])*sin(In[8]) + 0.0018798955148*cos(In[18])*cos(In[8]);
    Out[1+165] = 0.0014099216361*sin(In[19])*sin(In[8]) + 0.0014099216361*cos(In[19])*cos(In[8]);
    Out[1+166] = 0.0009399477574*sin(In[20])*sin(In[8]) + 0.0009399477574*cos(In[20])*cos(In[8]);
    Out[1+167] = 0.0004699738787*sin(In[21])*sin(In[8]) + 0.0004699738787*cos(In[21])*cos(In[8]);
    Out[1+168] = 0.103729379*cos(In[9]);
    Out[1+169] = 0.0061096604231*sin(In[2])*sin(In[9]) + 0.0061096604231*cos(In[2])*cos(In[9]);
    Out[1+170] = 0.0061096604231*sin(In[3])*sin(In[9]) + 0.0061096604231*cos(In[3])*cos(In[9]);
    Out[1+171] = 0.0061096604231*sin(In[4])*sin(In[9]) + 0.0061096604231*cos(In[4])*cos(In[9]);
    Out[1+172] = 0.0061096604231*sin(In[5])*sin(In[9]) + 0.0061096604231*cos(In[5])*cos(In[9]);
    Out[1+173] = 0.0061096604231*sin(In[6])*sin(In[9]) + 0.0061096604231*cos(In[6])*cos(In[9]);
    Out[1+174] = 0.0061096604231*sin(In[7])*sin(In[9]) + 0.0061096604231*cos(In[7])*cos(In[9]);
    Out[1+175] = 0.0061096604231*sin(In[8])*sin(In[9]) + 0.0061096604231*cos(In[8])*cos(In[9]);
    Out[1+176] = 0.0061096604231;
    Out[1+177] = 0.0056396865444*sin(In[10])*sin(In[9]) + 0.0056396865444*cos(In[10])*cos(In[9]);
    Out[1+178] = 0.0051697126657*sin(In[11])*sin(In[9]) + 0.0051697126657*cos(In[11])*cos(In[9]);
    Out[1+179] = 0.004699738787*sin(In[12])*sin(In[9]) + 0.004699738787*cos(In[12])*cos(In[9]);
    Out[1+180] = 0.0042297649083*sin(In[13])*sin(In[9]) + 0.0042297649083*cos(In[13])*cos(In[9]);
    Out[1+181] = 0.0037597910296*sin(In[14])*sin(In[9]) + 0.0037597910296*cos(In[14])*cos(In[9]);
    Out[1+182] = 0.0032898171509*sin(In[15])*sin(In[9]) + 0.0032898171509*cos(In[15])*cos(In[9]);
    Out[1+183] = 0.0028198432722*sin(In[16])*sin(In[9]) + 0.0028198432722*cos(In[16])*cos(In[9]);
    Out[1+184] = 0.0023498693935*sin(In[17])*sin(In[9]) + 0.0023498693935*cos(In[17])*cos(In[9]);
    Out[1+185] = 0.0018798955148*sin(In[18])*sin(In[9]) + 0.0018798955148*cos(In[18])*cos(In[9]);
    Out[1+186] = 0.0014099216361*sin(In[19])*sin(In[9]) + 0.0014099216361*cos(In[19])*cos(In[9]);
    Out[1+187] = 0.0009399477574*sin(In[20])*sin(In[9]) + 0.0009399477574*cos(In[20])*cos(In[9]);
    Out[1+188] = 0.0004699738787*sin(In[21])*sin(In[9]) + 0.0004699738787*cos(In[21])*cos(In[9]);
    Out[1+189] = 0.095750196*cos(In[10]);
    Out[1+190] = 0.0056396865444*sin(In[10])*sin(In[2]) + 0.0056396865444*cos(In[10])*cos(In[2]);
    Out[1+191] = 0.0056396865444*sin(In[10])*sin(In[3]) + 0.0056396865444*cos(In[10])*cos(In[3]);
    Out[1+192] = 0.0056396865444*sin(In[10])*sin(In[4]) + 0.0056396865444*cos(In[10])*cos(In[4]);
    Out[1+193] = 0.0056396865444*sin(In[10])*sin(In[5]) + 0.0056396865444*cos(In[10])*cos(In[5]);
    Out[1+194] = 0.0056396865444*sin(In[10])*sin(In[6]) + 0.0056396865444*cos(In[10])*cos(In[6]);
    Out[1+195] = 0.0056396865444*sin(In[10])*sin(In[7]) + 0.0056396865444*cos(In[10])*cos(In[7]);
    Out[1+196] = 0.0056396865444*sin(In[10])*sin(In[8]) + 0.0056396865444*cos(In[10])*cos(In[8]);
    Out[1+197] = 0.0056396865444*sin(In[10])*sin(In[9]) + 0.0056396865444*cos(In[10])*cos(In[9]);
    Out[1+198] = 0.0056396865444;
    Out[1+199] = 0.0051697126657*sin(In[10])*sin(In[11]) + 0.0051697126657*cos(In[10])*cos(In[11]);
    Out[1+200] = 0.004699738787*sin(In[10])*sin(In[12]) + 0.004699738787*cos(In[10])*cos(In[12]);
    Out[1+201] = 0.0042297649083*sin(In[10])*sin(In[13]) + 0.0042297649083*cos(In[10])*cos(In[13]);
    Out[1+202] = 0.0037597910296*sin(In[10])*sin(In[14]) + 0.0037597910296*cos(In[10])*cos(In[14]);
    Out[1+203] = 0.0032898171509*sin(In[10])*sin(In[15]) + 0.0032898171509*cos(In[10])*cos(In[15]);
    Out[1+204] = 0.0028198432722*sin(In[10])*sin(In[16]) + 0.0028198432722*cos(In[10])*cos(In[16]);
    Out[1+205] = 0.0023498693935*sin(In[10])*sin(In[17]) + 0.0023498693935*cos(In[10])*cos(In[17]);
    Out[1+206] = 0.0018798955148*sin(In[10])*sin(In[18]) + 0.0018798955148*cos(In[10])*cos(In[18]);
    Out[1+207] = 0.0014099216361*sin(In[10])*sin(In[19]) + 0.0014099216361*cos(In[10])*cos(In[19]);
    Out[1+208] = 0.0009399477574*sin(In[10])*sin(In[20]) + 0.0009399477574*cos(In[10])*cos(In[20]);
    Out[1+209] = 0.0004699738787*sin(In[10])*sin(In[21]) + 0.0004699738787*cos(In[10])*cos(In[21]);
    Out[1+210] = 0.087771013*cos(In[11]);
    Out[1+211] = 0.0051697126657*sin(In[11])*sin(In[2]) + 0.0051697126657*cos(In[11])*cos(In[2]);
    Out[1+212] = 0.0051697126657*sin(In[11])*sin(In[3]) + 0.0051697126657*cos(In[11])*cos(In[3]);
    Out[1+213] = 0.0051697126657*sin(In[11])*sin(In[4]) + 0.0051697126657*cos(In[11])*cos(In[4]);
    Out[1+214] = 0.0051697126657*sin(In[11])*sin(In[5]) + 0.0051697126657*cos(In[11])*cos(In[5]);
    Out[1+215] = 0.0051697126657*sin(In[11])*sin(In[6]) + 0.0051697126657*cos(In[11])*cos(In[6]);
    Out[1+216] = 0.0051697126657*sin(In[11])*sin(In[7]) + 0.0051697126657*cos(In[11])*cos(In[7]);
    Out[1+217] = 0.0051697126657*sin(In[11])*sin(In[8]) + 0.0051697126657*cos(In[11])*cos(In[8]);
    Out[1+218] = 0.0051697126657*sin(In[11])*sin(In[9]) + 0.0051697126657*cos(In[11])*cos(In[9]);
    Out[1+219] = 0.0051697126657*sin(In[10])*sin(In[11]) + 0.0051697126657*cos(In[10])*cos(In[11]);
    Out[1+220] = 0.0051697126657;
    Out[1+221] = 0.004699738787*sin(In[11])*sin(In[12]) + 0.004699738787*cos(In[11])*cos(In[12]);
    Out[1+222] = 0.0042297649083*sin(In[11])*sin(In[13]) + 0.0042297649083*cos(In[11])*cos(In[13]);
    Out[1+223] = 0.0037597910296*sin(In[11])*sin(In[14]) + 0.0037597910296*cos(In[11])*cos(In[14]);
    Out[1+224] = 0.0032898171509*sin(In[11])*sin(In[15]) + 0.0032898171509*cos(In[11])*cos(In[15]);
    Out[1+225] = 0.0028198432722*sin(In[11])*sin(In[16]) + 0.0028198432722*cos(In[11])*cos(In[16]);
    Out[1+226] = 0.0023498693935*sin(In[11])*sin(In[17]) + 0.0023498693935*cos(In[11])*cos(In[17]);
    Out[1+227] = 0.0018798955148*sin(In[11])*sin(In[18]) + 0.0018798955148*cos(In[11])*cos(In[18]);
    Out[1+228] = 0.0014099216361*sin(In[11])*sin(In[19]) + 0.0014099216361*cos(In[11])*cos(In[19]);
    Out[1+229] = 0.0009399477574*sin(In[11])*sin(In[20]) + 0.0009399477574*cos(In[11])*cos(In[20]);
    Out[1+230] = 0.0004699738787*sin(In[11])*sin(In[21]) + 0.0004699738787*cos(In[11])*cos(In[21]);
    Out[1+231] = 0.07979183*cos(In[12]);
    Out[1+232] = 0.004699738787*sin(In[12])*sin(In[2]) + 0.004699738787*cos(In[12])*cos(In[2]);
    Out[1+233] = 0.004699738787*sin(In[12])*sin(In[3]) + 0.004699738787*cos(In[12])*cos(In[3]);
    Out[1+234] = 0.004699738787*sin(In[12])*sin(In[4]) + 0.004699738787*cos(In[12])*cos(In[4]);
    Out[1+235] = 0.004699738787*sin(In[12])*sin(In[5]) + 0.004699738787*cos(In[12])*cos(In[5]);
    Out[1+236] = 0.004699738787*sin(In[12])*sin(In[6]) + 0.004699738787*cos(In[12])*cos(In[6]);
    Out[1+237] = 0.004699738787*sin(In[12])*sin(In[7]) + 0.004699738787*cos(In[12])*cos(In[7]);
    Out[1+238] = 0.004699738787*sin(In[12])*sin(In[8]) + 0.004699738787*cos(In[12])*cos(In[8]);
    Out[1+239] = 0.004699738787*sin(In[12])*sin(In[9]) + 0.004699738787*cos(In[12])*cos(In[9]);
    Out[1+240] = 0.004699738787*sin(In[10])*sin(In[12]) + 0.004699738787*cos(In[10])*cos(In[12]);
    Out[1+241] = 0.004699738787*sin(In[11])*sin(In[12]) + 0.004699738787*cos(In[11])*cos(In[12]);
    Out[1+242] = 0.004699738787;
    Out[1+243] = 0.0042297649083*sin(In[12])*sin(In[13]) + 0.0042297649083*cos(In[12])*cos(In[13]);
    Out[1+244] = 0.0037597910296*sin(In[12])*sin(In[14]) + 0.0037597910296*cos(In[12])*cos(In[14]);
    Out[1+245] = 0.0032898171509*sin(In[12])*sin(In[15]) + 0.0032898171509*cos(In[12])*cos(In[15]);
    Out[1+246] = 0.0028198432722*sin(In[12])*sin(In[16]) + 0.0028198432722*cos(In[12])*cos(In[16]);
    Out[1+247] = 0.0023498693935*sin(In[12])*sin(In[17]) + 0.0023498693935*cos(In[12])*cos(In[17]);
    Out[1+248] = 0.0018798955148*sin(In[12])*sin(In[18]) + 0.0018798955148*cos(In[12])*cos(In[18]);
    Out[1+249] = 0.0014099216361*sin(In[12])*sin(In[19]) + 0.0014099216361*cos(In[12])*cos(In[19]);
    Out[1+250] = 0.0009399477574*sin(In[12])*sin(In[20]) + 0.0009399477574*cos(In[12])*cos(In[20]);
    Out[1+251] = 0.0004699738787*sin(In[12])*sin(In[21]) + 0.0004699738787*cos(In[12])*cos(In[21]);
    Out[1+252] = 0.071812647*cos(In[13]);
    Out[1+253] = 0.0042297649083*sin(In[13])*sin(In[2]) + 0.0042297649083*cos(In[13])*cos(In[2]);
    Out[1+254] = 0.0042297649083*sin(In[13])*sin(In[3]) + 0.0042297649083*cos(In[13])*cos(In[3]);
    Out[1+255] = 0.0042297649083*sin(In[13])*sin(In[4]) + 0.0042297649083*cos(In[13])*cos(In[4]);
    Out[1+256] = 0.0042297649083*sin(In[13])*sin(In[5]) + 0.0042297649083*cos(In[13])*cos(In[5]);
    Out[1+257] = 0.0042297649083*sin(In[13])*sin(In[6]) + 0.0042297649083*cos(In[13])*cos(In[6]);
    Out[1+258] = 0.0042297649083*sin(In[13])*sin(In[7]) + 0.0042297649083*cos(In[13])*cos(In[7]);
    Out[1+259] = 0.0042297649083*sin(In[13])*sin(In[8]) + 0.0042297649083*cos(In[13])*cos(In[8]);
    Out[1+260] = 0.0042297649083*sin(In[13])*sin(In[9]) + 0.0042297649083*cos(In[13])*cos(In[9]);
    Out[1+261] = 0.0042297649083*sin(In[10])*sin(In[13]) + 0.0042297649083*cos(In[10])*cos(In[13]);
    Out[1+262] = 0.0042297649083*sin(In[11])*sin(In[13]) + 0.0042297649083*cos(In[11])*cos(In[13]);
    Out[1+263] = 0.0042297649083*sin(In[12])*sin(In[13]) + 0.0042297649083*cos(In[12])*cos(In[13]);
    Out[1+264] = 0.0042297649083;
    Out[1+265] = 0.0037597910296*sin(In[13])*sin(In[14]) + 0.0037597910296*cos(In[13])*cos(In[14]);
    Out[1+266] = 0.0032898171509*sin(In[13])*sin(In[15]) + 0.0032898171509*cos(In[13])*cos(In[15]);
    Out[1+267] = 0.0028198432722*sin(In[13])*sin(In[16]) + 0.0028198432722*cos(In[13])*cos(In[16]);
    Out[1+268] = 0.0023498693935*sin(In[13])*sin(In[17]) + 0.0023498693935*cos(In[13])*cos(In[17]);
    Out[1+269] = 0.0018798955148*sin(In[13])*sin(In[18]) + 0.0018798955148*cos(In[13])*cos(In[18]);
    Out[1+270] = 0.0014099216361*sin(In[13])*sin(In[19]) + 0.0014099216361*cos(In[13])*cos(In[19]);
    Out[1+271] = 0.0009399477574*sin(In[13])*sin(In[20]) + 0.0009399477574*cos(In[13])*cos(In[20]);
    Out[1+272] = 0.0004699738787*sin(In[13])*sin(In[21]) + 0.0004699738787*cos(In[13])*cos(In[21]);
    Out[1+273] = 0.063833464*cos(In[14]);
    Out[1+274] = 0.0037597910296*sin(In[14])*sin(In[2]) + 0.0037597910296*cos(In[14])*cos(In[2]);
    Out[1+275] = 0.0037597910296*sin(In[14])*sin(In[3]) + 0.0037597910296*cos(In[14])*cos(In[3]);
    Out[1+276] = 0.0037597910296*sin(In[14])*sin(In[4]) + 0.0037597910296*cos(In[14])*cos(In[4]);
    Out[1+277] = 0.0037597910296*sin(In[14])*sin(In[5]) + 0.0037597910296*cos(In[14])*cos(In[5]);
    Out[1+278] = 0.0037597910296*sin(In[14])*sin(In[6]) + 0.0037597910296*cos(In[14])*cos(In[6]);
    Out[1+279] = 0.0037597910296*sin(In[14])*sin(In[7]) + 0.0037597910296*cos(In[14])*cos(In[7]);
    Out[1+280] = 0.0037597910296*sin(In[14])*sin(In[8]) + 0.0037597910296*cos(In[14])*cos(In[8]);
    Out[1+281] = 0.0037597910296*sin(In[14])*sin(In[9]) + 0.0037597910296*cos(In[14])*cos(In[9]);
    Out[1+282] = 0.0037597910296*sin(In[10])*sin(In[14]) + 0.0037597910296*cos(In[10])*cos(In[14]);
    Out[1+283] = 0.0037597910296*sin(In[11])*sin(In[14]) + 0.0037597910296*cos(In[11])*cos(In[14]);
    Out[1+284] = 0.0037597910296*sin(In[12])*sin(In[14]) + 0.0037597910296*cos(In[12])*cos(In[14]);
    Out[1+285] = 0.0037597910296*sin(In[13])*sin(In[14]) + 0.0037597910296*cos(In[13])*cos(In[14]);
    Out[1+286] = 0.0037597910296;
    Out[1+287] = 0.0032898171509*sin(In[14])*sin(In[15]) + 0.0032898171509*cos(In[14])*cos(In[15]);
    Out[1+288] = 0.0028198432722*sin(In[14])*sin(In[16]) + 0.0028198432722*cos(In[14])*cos(In[16]);
    Out[1+289] = 0.0023498693935*sin(In[14])*sin(In[17]) + 0.0023498693935*cos(In[14])*cos(In[17]);
    Out[1+290] = 0.0018798955148*sin(In[14])*sin(In[18]) + 0.0018798955148*cos(In[14])*cos(In[18]);
    Out[1+291] = 0.0014099216361*sin(In[14])*sin(In[19]) + 0.0014099216361*cos(In[14])*cos(In[19]);
    Out[1+292] = 0.0009399477574*sin(In[14])*sin(In[20]) + 0.0009399477574*cos(In[14])*cos(In[20]);
    Out[1+293] = 0.0004699738787*sin(In[14])*sin(In[21]) + 0.0004699738787*cos(In[14])*cos(In[21]);
    Out[1+294] = 0.055854281*cos(In[15]);
    Out[1+295] = 0.0032898171509*sin(In[15])*sin(In[2]) + 0.0032898171509*cos(In[15])*cos(In[2]);
    Out[1+296] = 0.0032898171509*sin(In[15])*sin(In[3]) + 0.0032898171509*cos(In[15])*cos(In[3]);
    Out[1+297] = 0.0032898171509*sin(In[15])*sin(In[4]) + 0.0032898171509*cos(In[15])*cos(In[4]);
    Out[1+298] = 0.0032898171509*sin(In[15])*sin(In[5]) + 0.0032898171509*cos(In[15])*cos(In[5]);
    Out[1+299] = 0.0032898171509*sin(In[15])*sin(In[6]) + 0.0032898171509*cos(In[15])*cos(In[6]);
    Out[1+300] = 0.0032898171509*sin(In[15])*sin(In[7]) + 0.0032898171509*cos(In[15])*cos(In[7]);
    Out[1+301] = 0.0032898171509*sin(In[15])*sin(In[8]) + 0.0032898171509*cos(In[15])*cos(In[8]);
    Out[1+302] = 0.0032898171509*sin(In[15])*sin(In[9]) + 0.0032898171509*cos(In[15])*cos(In[9]);
    Out[1+303] = 0.0032898171509*sin(In[10])*sin(In[15]) + 0.0032898171509*cos(In[10])*cos(In[15]);
    Out[1+304] = 0.0032898171509*sin(In[11])*sin(In[15]) + 0.0032898171509*cos(In[11])*cos(In[15]);
    Out[1+305] = 0.0032898171509*sin(In[12])*sin(In[15]) + 0.0032898171509*cos(In[12])*cos(In[15]);
    Out[1+306] = 0.0032898171509*sin(In[13])*sin(In[15]) + 0.0032898171509*cos(In[13])*cos(In[15]);
    Out[1+307] = 0.0032898171509*sin(In[14])*sin(In[15]) + 0.0032898171509*cos(In[14])*cos(In[15]);
    Out[1+308] = 0.0032898171509;
    Out[1+309] = 0.0028198432722*sin(In[15])*sin(In[16]) + 0.0028198432722*cos(In[15])*cos(In[16]);
    Out[1+310] = 0.0023498693935*sin(In[15])*sin(In[17]) + 0.0023498693935*cos(In[15])*cos(In[17]);
    Out[1+311] = 0.0018798955148*sin(In[15])*sin(In[18]) + 0.0018798955148*cos(In[15])*cos(In[18]);
    Out[1+312] = 0.0014099216361*sin(In[15])*sin(In[19]) + 0.0014099216361*cos(In[15])*cos(In[19]);
    Out[1+313] = 0.0009399477574*sin(In[15])*sin(In[20]) + 0.0009399477574*cos(In[15])*cos(In[20]);
    Out[1+314] = 0.0004699738787*sin(In[15])*sin(In[21]) + 0.0004699738787*cos(In[15])*cos(In[21]);
    Out[1+315] = 0.047875098*cos(In[16]);
    Out[1+316] = 0.0028198432722*sin(In[16])*sin(In[2]) + 0.0028198432722*cos(In[16])*cos(In[2]);
    Out[1+317] = 0.0028198432722*sin(In[16])*sin(In[3]) + 0.0028198432722*cos(In[16])*cos(In[3]);
    Out[1+318] = 0.0028198432722*sin(In[16])*sin(In[4]) + 0.0028198432722*cos(In[16])*cos(In[4]);
    Out[1+319] = 0.0028198432722*sin(In[16])*sin(In[5]) + 0.0028198432722*cos(In[16])*cos(In[5]);
    Out[1+320] = 0.0028198432722*sin(In[16])*sin(In[6]) + 0.0028198432722*cos(In[16])*cos(In[6]);
    Out[1+321] = 0.0028198432722*sin(In[16])*sin(In[7]) + 0.0028198432722*cos(In[16])*cos(In[7]);
    Out[1+322] = 0.0028198432722*sin(In[16])*sin(In[8]) + 0.0028198432722*cos(In[16])*cos(In[8]);
    Out[1+323] = 0.0028198432722*sin(In[16])*sin(In[9]) + 0.0028198432722*cos(In[16])*cos(In[9]);
    Out[1+324] = 0.0028198432722*sin(In[10])*sin(In[16]) + 0.0028198432722*cos(In[10])*cos(In[16]);
    Out[1+325] = 0.0028198432722*sin(In[11])*sin(In[16]) + 0.0028198432722*cos(In[11])*cos(In[16]);
    Out[1+326] = 0.0028198432722*sin(In[12])*sin(In[16]) + 0.0028198432722*cos(In[12])*cos(In[16]);
    Out[1+327] = 0.0028198432722*sin(In[13])*sin(In[16]) + 0.0028198432722*cos(In[13])*cos(In[16]);
    Out[1+328] = 0.0028198432722*sin(In[14])*sin(In[16]) + 0.0028198432722*cos(In[14])*cos(In[16]);
    Out[1+329] = 0.0028198432722*sin(In[15])*sin(In[16]) + 0.0028198432722*cos(In[15])*cos(In[16]);
    Out[1+330] = 0.0028198432722;
    Out[1+331] = 0.0023498693935*sin(In[16])*sin(In[17]) + 0.0023498693935*cos(In[16])*cos(In[17]);
    Out[1+332] = 0.0018798955148*sin(In[16])*sin(In[18]) + 0.0018798955148*cos(In[16])*cos(In[18]);
    Out[1+333] = 0.0014099216361*sin(In[16])*sin(In[19]) + 0.0014099216361*cos(In[16])*cos(In[19]);
    Out[1+334] = 0.0009399477574*sin(In[16])*sin(In[20]) + 0.0009399477574*cos(In[16])*cos(In[20]);
    Out[1+335] = 0.0004699738787*sin(In[16])*sin(In[21]) + 0.0004699738787*cos(In[16])*cos(In[21]);
    Out[1+336] = 0.039895915*cos(In[17]);
    Out[1+337] = 0.0023498693935*sin(In[17])*sin(In[2]) + 0.0023498693935*cos(In[17])*cos(In[2]);
    Out[1+338] = 0.0023498693935*sin(In[17])*sin(In[3]) + 0.0023498693935*cos(In[17])*cos(In[3]);
    Out[1+339] = 0.0023498693935*sin(In[17])*sin(In[4]) + 0.0023498693935*cos(In[17])*cos(In[4]);
    Out[1+340] = 0.0023498693935*sin(In[17])*sin(In[5]) + 0.0023498693935*cos(In[17])*cos(In[5]);
    Out[1+341] = 0.0023498693935*sin(In[17])*sin(In[6]) + 0.0023498693935*cos(In[17])*cos(In[6]);
    Out[1+342] = 0.0023498693935*sin(In[17])*sin(In[7]) + 0.0023498693935*cos(In[17])*cos(In[7]);
    Out[1+343] = 0.0023498693935*sin(In[17])*sin(In[8]) + 0.0023498693935*cos(In[17])*cos(In[8]);
    Out[1+344] = 0.0023498693935*sin(In[17])*sin(In[9]) + 0.0023498693935*cos(In[17])*cos(In[9]);
    Out[1+345] = 0.0023498693935*sin(In[10])*sin(In[17]) + 0.0023498693935*cos(In[10])*cos(In[17]);
    Out[1+346] = 0.0023498693935*sin(In[11])*sin(In[17]) + 0.0023498693935*cos(In[11])*cos(In[17]);
    Out[1+347] = 0.0023498693935*sin(In[12])*sin(In[17]) + 0.0023498693935*cos(In[12])*cos(In[17]);
    Out[1+348] = 0.0023498693935*sin(In[13])*sin(In[17]) + 0.0023498693935*cos(In[13])*cos(In[17]);
    Out[1+349] = 0.0023498693935*sin(In[14])*sin(In[17]) + 0.0023498693935*cos(In[14])*cos(In[17]);
    Out[1+350] = 0.0023498693935*sin(In[15])*sin(In[17]) + 0.0023498693935*cos(In[15])*cos(In[17]);
    Out[1+351] = 0.0023498693935*sin(In[16])*sin(In[17]) + 0.0023498693935*cos(In[16])*cos(In[17]);
    Out[1+352] = 0.0023498693935;
    Out[1+353] = 0.0018798955148*sin(In[17])*sin(In[18]) + 0.0018798955148*cos(In[17])*cos(In[18]);
    Out[1+354] = 0.0014099216361*sin(In[17])*sin(In[19]) + 0.0014099216361*cos(In[17])*cos(In[19]);
    Out[1+355] = 0.0009399477574*sin(In[17])*sin(In[20]) + 0.0009399477574*cos(In[17])*cos(In[20]);
    Out[1+356] = 0.0004699738787*sin(In[17])*sin(In[21]) + 0.0004699738787*cos(In[17])*cos(In[21]);
    Out[1+357] = 0.031916732*cos(In[18]);
    Out[1+358] = 0.0018798955148*sin(In[18])*sin(In[2]) + 0.0018798955148*cos(In[18])*cos(In[2]);
    Out[1+359] = 0.0018798955148*sin(In[18])*sin(In[3]) + 0.0018798955148*cos(In[18])*cos(In[3]);
    Out[1+360] = 0.0018798955148*sin(In[18])*sin(In[4]) + 0.0018798955148*cos(In[18])*cos(In[4]);
    Out[1+361] = 0.0018798955148*sin(In[18])*sin(In[5]) + 0.0018798955148*cos(In[18])*cos(In[5]);
    Out[1+362] = 0.0018798955148*sin(In[18])*sin(In[6]) + 0.0018798955148*cos(In[18])*cos(In[6]);
    Out[1+363] = 0.0018798955148*sin(In[18])*sin(In[7]) + 0.0018798955148*cos(In[18])*cos(In[7]);
    Out[1+364] = 0.0018798955148*sin(In[18])*sin(In[8]) + 0.0018798955148*cos(In[18])*cos(In[8]);
    Out[1+365] = 0.0018798955148*sin(In[18])*sin(In[9]) + 0.0018798955148*cos(In[18])*cos(In[9]);
    Out[1+366] = 0.0018798955148*sin(In[10])*sin(In[18]) + 0.0018798955148*cos(In[10])*cos(In[18]);
    Out[1+367] = 0.0018798955148*sin(In[11])*sin(In[18]) + 0.0018798955148*cos(In[11])*cos(In[18]);
    Out[1+368] = 0.0018798955148*sin(In[12])*sin(In[18]) + 0.0018798955148*cos(In[12])*cos(In[18]);
    Out[1+369] = 0.0018798955148*sin(In[13])*sin(In[18]) + 0.0018798955148*cos(In[13])*cos(In[18]);
    Out[1+370] = 0.0018798955148*sin(In[14])*sin(In[18]) + 0.0018798955148*cos(In[14])*cos(In[18]);
    Out[1+371] = 0.0018798955148*sin(In[15])*sin(In[18]) + 0.0018798955148*cos(In[15])*cos(In[18]);
    Out[1+372] = 0.0018798955148*sin(In[16])*sin(In[18]) + 0.0018798955148*cos(In[16])*cos(In[18]);
    Out[1+373] = 0.0018798955148*sin(In[17])*sin(In[18]) + 0.0018798955148*cos(In[17])*cos(In[18]);
    Out[1+374] = 0.0018798955148;
    Out[1+375] = 0.0014099216361*sin(In[18])*sin(In[19]) + 0.0014099216361*cos(In[18])*cos(In[19]);
    Out[1+376] = 0.0009399477574*sin(In[18])*sin(In[20]) + 0.0009399477574*cos(In[18])*cos(In[20]);
    Out[1+377] = 0.0004699738787*sin(In[18])*sin(In[21]) + 0.0004699738787*cos(In[18])*cos(In[21]);
    Out[1+378] = 0.023937549*cos(In[19]);
    Out[1+379] = 0.0014099216361*sin(In[19])*sin(In[2]) + 0.0014099216361*cos(In[19])*cos(In[2]);
    Out[1+380] = 0.0014099216361*sin(In[19])*sin(In[3]) + 0.0014099216361*cos(In[19])*cos(In[3]);
    Out[1+381] = 0.0014099216361*sin(In[19])*sin(In[4]) + 0.0014099216361*cos(In[19])*cos(In[4]);
    Out[1+382] = 0.0014099216361*sin(In[19])*sin(In[5]) + 0.0014099216361*cos(In[19])*cos(In[5]);
    Out[1+383] = 0.0014099216361*sin(In[19])*sin(In[6]) + 0.0014099216361*cos(In[19])*cos(In[6]);
    Out[1+384] = 0.0014099216361*sin(In[19])*sin(In[7]) + 0.0014099216361*cos(In[19])*cos(In[7]);
    Out[1+385] = 0.0014099216361*sin(In[19])*sin(In[8]) + 0.0014099216361*cos(In[19])*cos(In[8]);
    Out[1+386] = 0.0014099216361*sin(In[19])*sin(In[9]) + 0.0014099216361*cos(In[19])*cos(In[9]);
    Out[1+387] = 0.0014099216361*sin(In[10])*sin(In[19]) + 0.0014099216361*cos(In[10])*cos(In[19]);
    Out[1+388] = 0.0014099216361*sin(In[11])*sin(In[19]) + 0.0014099216361*cos(In[11])*cos(In[19]);
    Out[1+389] = 0.0014099216361*sin(In[12])*sin(In[19]) + 0.0014099216361*cos(In[12])*cos(In[19]);
    Out[1+390] = 0.0014099216361*sin(In[13])*sin(In[19]) + 0.0014099216361*cos(In[13])*cos(In[19]);
    Out[1+391] = 0.0014099216361*sin(In[14])*sin(In[19]) + 0.0014099216361*cos(In[14])*cos(In[19]);
    Out[1+392] = 0.0014099216361*sin(In[15])*sin(In[19]) + 0.0014099216361*cos(In[15])*cos(In[19]);
    Out[1+393] = 0.0014099216361*sin(In[16])*sin(In[19]) + 0.0014099216361*cos(In[16])*cos(In[19]);
    Out[1+394] = 0.0014099216361*sin(In[17])*sin(In[19]) + 0.0014099216361*cos(In[17])*cos(In[19]);
    Out[1+395] = 0.0014099216361*sin(In[18])*sin(In[19]) + 0.0014099216361*cos(In[18])*cos(In[19]);
    Out[1+396] = 0.0014099216361;
    Out[1+397] = 0.0009399477574*sin(In[19])*sin(In[20]) + 0.0009399477574*cos(In[19])*cos(In[20]);
    Out[1+398] = 0.0004699738787*sin(In[19])*sin(In[21]) + 0.0004699738787*cos(In[19])*cos(In[21]);
    Out[1+399] = 0.015958366*cos(In[20]);
    Out[1+400] = 0.0009399477574*sin(In[20])*sin(In[2]) + 0.0009399477574*cos(In[20])*cos(In[2]);
    Out[1+401] = 0.0009399477574*sin(In[20])*sin(In[3]) + 0.0009399477574*cos(In[20])*cos(In[3]);
    Out[1+402] = 0.0009399477574*sin(In[20])*sin(In[4]) + 0.0009399477574*cos(In[20])*cos(In[4]);
    Out[1+403] = 0.0009399477574*sin(In[20])*sin(In[5]) + 0.0009399477574*cos(In[20])*cos(In[5]);
    Out[1+404] = 0.0009399477574*sin(In[20])*sin(In[6]) + 0.0009399477574*cos(In[20])*cos(In[6]);
    Out[1+405] = 0.0009399477574*sin(In[20])*sin(In[7]) + 0.0009399477574*cos(In[20])*cos(In[7]);
    Out[1+406] = 0.0009399477574*sin(In[20])*sin(In[8]) + 0.0009399477574*cos(In[20])*cos(In[8]);
    Out[1+407] = 0.0009399477574*sin(In[20])*sin(In[9]) + 0.0009399477574*cos(In[20])*cos(In[9]);
    Out[1+408] = 0.0009399477574*sin(In[10])*sin(In[20]) + 0.0009399477574*cos(In[10])*cos(In[20]);
    Out[1+409] = 0.0009399477574*sin(In[11])*sin(In[20]) + 0.0009399477574*cos(In[11])*cos(In[20]);
    Out[1+410] = 0.0009399477574*sin(In[12])*sin(In[20]) + 0.0009399477574*cos(In[12])*cos(In[20]);
    Out[1+411] = 0.0009399477574*sin(In[13])*sin(In[20]) + 0.0009399477574*cos(In[13])*cos(In[20]);
    Out[1+412] = 0.0009399477574*sin(In[14])*sin(In[20]) + 0.0009399477574*cos(In[14])*cos(In[20]);
    Out[1+413] = 0.0009399477574*sin(In[15])*sin(In[20]) + 0.0009399477574*cos(In[15])*cos(In[20]);
    Out[1+414] = 0.0009399477574*sin(In[16])*sin(In[20]) + 0.0009399477574*cos(In[16])*cos(In[20]);
    Out[1+415] = 0.0009399477574*sin(In[17])*sin(In[20]) + 0.0009399477574*cos(In[17])*cos(In[20]);
    Out[1+416] = 0.0009399477574*sin(In[18])*sin(In[20]) + 0.0009399477574*cos(In[18])*cos(In[20]);
    Out[1+417] = 0.0009399477574*sin(In[19])*sin(In[20]) + 0.0009399477574*cos(In[19])*cos(In[20]);
    Out[1+418] = 0.0009399477574;
    Out[1+419] = 0.0004699738787*sin(In[20])*sin(In[21]) + 0.0004699738787*cos(In[20])*cos(In[21]);
    Out[1+420] = 0.007979183*cos(In[21]);
    Out[1+421] = 0.0004699738787*sin(In[21])*sin(In[2]) + 0.0004699738787*cos(In[21])*cos(In[2]);
    Out[1+422] = 0.0004699738787*sin(In[21])*sin(In[3]) + 0.0004699738787*cos(In[21])*cos(In[3]);
    Out[1+423] = 0.0004699738787*sin(In[21])*sin(In[4]) + 0.0004699738787*cos(In[21])*cos(In[4]);
    Out[1+424] = 0.0004699738787*sin(In[21])*sin(In[5]) + 0.0004699738787*cos(In[21])*cos(In[5]);
    Out[1+425] = 0.0004699738787*sin(In[21])*sin(In[6]) + 0.0004699738787*cos(In[21])*cos(In[6]);
    Out[1+426] = 0.0004699738787*sin(In[21])*sin(In[7]) + 0.0004699738787*cos(In[21])*cos(In[7]);
    Out[1+427] = 0.0004699738787*sin(In[21])*sin(In[8]) + 0.0004699738787*cos(In[21])*cos(In[8]);
    Out[1+428] = 0.0004699738787*sin(In[21])*sin(In[9]) + 0.0004699738787*cos(In[21])*cos(In[9]);
    Out[1+429] = 0.0004699738787*sin(In[10])*sin(In[21]) + 0.0004699738787*cos(In[10])*cos(In[21]);
    Out[1+430] = 0.0004699738787*sin(In[11])*sin(In[21]) + 0.0004699738787*cos(In[11])*cos(In[21]);
    Out[1+431] = 0.0004699738787*sin(In[12])*sin(In[21]) + 0.0004699738787*cos(In[12])*cos(In[21]);
    Out[1+432] = 0.0004699738787*sin(In[13])*sin(In[21]) + 0.0004699738787*cos(In[13])*cos(In[21]);
    Out[1+433] = 0.0004699738787*sin(In[14])*sin(In[21]) + 0.0004699738787*cos(In[14])*cos(In[21]);
    Out[1+434] = 0.0004699738787*sin(In[15])*sin(In[21]) + 0.0004699738787*cos(In[15])*cos(In[21]);
    Out[1+435] = 0.0004699738787*sin(In[16])*sin(In[21]) + 0.0004699738787*cos(In[16])*cos(In[21]);
    Out[1+436] = 0.0004699738787*sin(In[17])*sin(In[21]) + 0.0004699738787*cos(In[17])*cos(In[21]);
    Out[1+437] = 0.0004699738787*sin(In[18])*sin(In[21]) + 0.0004699738787*cos(In[18])*cos(In[21]);
    Out[1+438] = 0.0004699738787*sin(In[19])*sin(In[21]) + 0.0004699738787*cos(In[19])*cos(In[21]);
    Out[1+439] = 0.0004699738787*sin(In[20])*sin(In[21]) + 0.0004699738787*cos(In[20])*cos(In[21]);
    Out[1+440] = 0.0004699738787;

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using FunctionOperators, Test

@testset "Show.jl" begin
    Op₁ = FunctionOperator{Float64}("Op₁", x -> x .^ 3, x -> cbrt.(x), (10,10), (10,10))
    @test sprint(show, Op₁) == "FunctionOperator{Float64}(Op₁, (10, 10), (10, 10))"
    @test sprint(show, "text/plain", Op₁, context = :module=>@__MODULE__) == "FunctionOperator with eltype Float64\n    Name: Op₁\n    Input dimensions: (10, 10)\n    Output dimensions: (10, 10)"
    @test sprint(show, Op₁ * Op₁) == "FunctionOperatorComposite{Float64}(Op₁ * Op₁, (10, 10), (10, 10), no plan)"
    combined = Op₁ * Op₁
    combined * ones(10, 10)
    @test sprint(show, combined) == "FunctionOperatorComposite{Float64}(Op₁ * Op₁, (10, 10), (10, 10), buffer0 .= Op₁.forw(buffer0 .= Op₁.forw(x)))"
    @test sprint(show, "text/plain", Op₁ * Op₁, context = :module=>@__MODULE__) == "FunctionOperatorComposite with eltype Float64\n    Name: Op₁ * Op₁\n    Input dimensions: (10, 10)\n    Output dimensions: (10, 10)\n    Plan: no plan"
end


text:
module ImageCoreTests

using ImageCore, Test, ReferenceTests

# If we've run the tests previously, there might be ambiguities from other packages
if :StatsBase ∉ map(x->Symbol(string(x)), values(Base.loaded_modules))
    @test isempty(detect_ambiguities(ImageCore, Base, Core))
end

using Documenter
DocMeta.setdocmeta!(ImageCore, :DocTestSetup, :(using ImageCore); recursive=true)
doctest(ImageCore, manual = false)

include("colorchannels.jl")
include("views.jl")
include("convert_reinterpret.jl")
include("traits.jl")
include("map.jl")
include("functions.jl")
include("show.jl")

# run these last
isCI = haskey(ENV, "CI") || get(ENV, "JULIA_PKGEVAL", false)
if Base.JLOptions().can_inline == 1 && !isCI
    @info "running benchmarks"
    include("benchmarks.jl")  # these fail if inlining is off
end

end


text:
module CellValuesOperationsTests

using Test
using TensorValues
using Gridap

using Gridap.CachedArrays

using ..CellValuesMocks
using ..MapsMocks

l = 10

# CellNumber Operations

ax = [1.2, VectorValue(1.3,2.0), VectorValue(1.3,2.0)]
bx = [1  , VectorValue(1.1,2.2), 1                   ]

for op in (:+,:-)
  @eval begin

    for a in ax
      v = TestIterCellValue(a,l)
      w = $op(v)
      o = [ $op(vi) for vi in v ]
      test_iter_cell_value(w,o)
    end

    for a in ax
      v = TestIndexCellValue(a,l)
      w = $op(v)
      o = [ $op(vi) for vi in v ]
      test_index_cell_value(w,o)
    end

    for (a,b) in zip(ax,bx)
      u = TestIterCellValue(a,l)
      v = TestIterCellValue(b,l)
      w = $op(u,v)
      o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
      test_iter_cell_value(w,o)
    end

    for (a,b) in zip(ax,bx)
      u = TestIndexCellValue(a,l)
      v = TestIndexCellValue(b,l)
      w = $op(u,v)
      o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
      test_index_cell_value(w,o)
    end

  end
end

a = TensorValue(1.0,0.0,1.0,2.0)

for op in (:inv,:det,:meas,:transpose)
  @eval begin

    v = TestIterCellValue(a,l)
    w = $op(v)
    o = [ $op(vi) for vi in v ]
    test_iter_cell_value(w,o)

    v = TestIndexCellValue(a,l)
    w = $op(v)
    o = [ $op(vi) for vi in v ]
    test_index_cell_value(w,o)

  end
end

a = VectorValue(1.0,0.0,1.0)
b = VectorValue(1.0,4.0,-1.0)

for op in (:inner,:outer)
  @eval begin

    u = TestIterCellValue(a,l)
    v = TestIterCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_iter_cell_value(w,o)

    u = TestIndexCellValue(a,l)
    v = TestIndexCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_index_cell_value(w,o)

  end
end

a = TensorValue(1.0,0.0,1.0,2.0)
b = VectorValue(1.0,4.0)

for op in (:*,:\)
  @eval begin

    u = TestIterCellValue(a,l)
    v = TestIterCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_iter_cell_value(w,o)

    u = TestIndexCellValue(a,l)
    v = TestIndexCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_index_cell_value(w,o)

  end
end

a = TensorValue(1.0,0.0,1.0,2.0)
b = 2

for op in (:/,)
  @eval begin

    u = TestIterCellValue(a,l)
    v = TestIterCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_iter_cell_value(w,o)

    u = TestIndexCellValue(a,l)
    v = TestIndexCellValue(b,l)
    w = $op(u,v)
    o = [ $op(ui,vi) for (ui,vi) in zip(u,v) ]
    test_index_cell_value(w,o)

  end
end

u = TestIterCellValue(a,l)
v = TestIterCellValue(a,l)
@test u == v
@test u ≈ v

# CellArray Operations

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)

ax = [rand(1,3,4), 1.0      , [v1,v2,v3], [v1,v2,v3], 1         ]
bx = [rand(2,3,1), rand(2,3), [v2,v3,v1], v1        , [v2,v3,v1]]
cx = [rand(1,3,4), [v1,v2,v3]]

for op in (:+,:-)
  @eval begin

    for a in cx

      v = TestIterCellValue(a,l)
      w = $op(v)
      o = [ CachedArray(broadcast($op,vi)) for vi in v ]
      test_iter_cell_array(w,o)

      v = TestIndexCellValue(a,l)
      w = $op(v)
      o = [ CachedArray(broadcast($op,vi)) for vi in v ]
      test_index_cell_array(w,o)

    end

    for (a,b) in zip(ax,bx)

      u = TestIterCellValue(a,l)
      v = TestIterCellValue(b,l)
      o = [ CachedArray(broadcast($op,ui,vi)) for (ui,vi) in zip(u,v) ]
      w = $op(u,v)
      test_iter_cell_array(w,o)

      u = TestIndexCellValue(a,l)
      v = TestIndexCellValue(b,l)
      o = [ CachedArray(broadcast($op,ui,vi)) for (ui,vi) in zip(u,v) ]
      w = $op(u,v)
      test_index_cell_array(w,o)

    end

  end
end

ai = TensorValue(1.0,0.0,1.0,2.0)
a = fill(ai,2,3)

for op in (:inv,:det,:meas,:transpose)
  @eval begin

  v = TestIterCellValue(a,l)
  w = $op(v)
  o = [ CachedArray(broadcast($op,vi)) for vi in v ]
  test_iter_cell_array(w,o)

  v = TestIndexCellValue(a,l)
  w = $op(v)
  o = [ CachedArray(broadcast($op,vi)) for vi in v ]
  test_index_cell_array(w,o)

  end
end

ai = VectorValue(1.0,0.0,1.0)
bi = VectorValue(1.0,4.0,-1.0)
a = fill(ai,2,3)
b = fill(bi,1,3)

for op in (:inner,:outer)
  @eval begin

    u = TestIterCellValue(a,l)
    v = TestIterCellValue(b,l)
    o = [ CachedArray(broadcast($op,ui,vi)) for (ui,vi) in zip(u,v) ]
    w = $op(u,v)
    test_iter_cell_array(w,o)

    u = TestIndexCellValue(a,l)
    v = TestIndexCellValue(b,l)
    o = [ CachedArray(broadcast($op,ui,vi)) for (ui,vi) in zip(u,v) ]
    w = $op(u,v)
    test_index_cell_array(w,o)

  end
end

# Eq

u = TestIterCellValue(a,l)
v = TestIterCellValue(a,l)
@test u == v
@test u ≈ v

# cellsum

D = 3
a = rand(2,3,4)
r = reshape(sum(a,dims=D),(2,3))
o = [CachedArray(r) for i in 1:l]

u = TestIterCellValue(a,l)
w = cellsum(u,dim=D)
test_iter_cell_array(w,o)

u = TestIndexCellValue(a,l)
w = cellsum(u,dim=D)
test_index_cell_array(w,o)

a = rand(4)
r = sum(a)
o = [r for i in 1:l]

u = TestIterCellValue(a,l)
w = cellsum(u,dim=1)
@test isa(w,CellNumber)
test_iter_cell_value(w,o)

# cellnewaxis

a = rand(2,3,4)
r = collect(reshape(a,(1,2,3,4)))
o = [CachedArray(r) for i in 1:l]

u = TestIterCellValue(a,l)
w = cellnewaxis(u,dim=1)
test_iter_cell_array(w,o)

u = TestIndexCellValue(a,l)
w = cellnewaxis(u,dim=1)
test_index_cell_array(w,o)

# cellmean

a = rand(2,3,4)
r = sum(a) / length(a)
o = [r for i in 1:l]

u = TestIterCellValue(a,l)
w = cellmean(u)
@test isa(w,CellNumber)
test_iter_cell_value(w,o)

# compress

a = [1,3,2]
u = TestIterCellValue(a,5)
u_data, u_ptrs = compress(u)
@test u_data == [1, 3, 2, 1, 3, 2, 1, 3, 2, 1, 3, 2, 1, 3, 2]
@test u_ptrs == [1, 4, 7, 10, 13, 16]

# CellMaps Operations

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)
r = evaluate(m,p)

for op in (:+,:-)
  @eval begin

    cm = TestIterCellValue(m,l)
    cp = TestIterCellValue(p,l)
    rm = [ CachedArray(broadcast($op,r)) for i in 1:l]
    cm2 = $op(cm)
    test_iter_cell_map(cm2,cp,rm)

    cm = TestIndexCellValue(m,l)
    cp = TestIndexCellValue(p,l)
    rm = [ CachedArray(broadcast($op,r)) for i in 1:l]
    cm2 = $op(cm)
    ca = evaluate(cm2,cp)
    test_index_cell_array(ca,rm)

    z = 2.0
    cm = TestIndexCellValue(m,l)
    cp = TestIndexCellValue(p,l)
    rm = [ CachedArray(broadcast($op,z,r)) for i in 1:l]
    cm2 = $op(z,cm)
    ca = evaluate(cm2,cp)
    test_index_cell_array(ca,rm)

    z = 2.0
    cm = TestIndexCellValue(m,l)
    cp = TestIndexCellValue(p,l)
    rm = [ CachedArray(broadcast($op,r,z)) for i in 1:l]
    cm2 = $op(cm,z)
    ca = evaluate(cm2,cp)
    test_index_cell_array(ca,rm)

  end
end

for op in (:inner,:outer)
  @eval begin

    cm = TestIterCellValue(m,l)
    cp = TestIterCellValue(p,l)
    rm = [ CachedArray(broadcast($op,r,r)) for i in 1:l]
    cm2 = $op(cm,cm)
    test_iter_cell_map(cm2,cp,rm)

    cm = TestIndexCellValue(m,l)
    cp = TestIndexCellValue(p,l)
    rm = [ CachedArray(broadcast($op,r,r)) for i in 1:l]
    cm2 = $op(cm,cm)
    ca = evaluate(cm2,cp)
    test_index_cell_array(ca,rm)

  end
end

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rc = collect(reshape(r,(1,3)))
rm = [ CachedArray(rc) for i in 1:l]
cm2 = cellnewaxis(cm,dim=1)
test_iter_cell_map(cm2,cp,rm)

end # module


text:
function Bias(arg0::Vector{JString}, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::Vector{jdouble}, arg4::Vector{jdouble})
    return Bias((Vector{JString}, Vector{jdouble}, Vector{jdouble}, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4)
end

function get_parameters_drivers(obj::Bias)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::Bias, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
# Set up problem using TrajectoryOptimization.jl and RobotZoo.jl
using TrajectoryOptimization
using Altro
import RobotZoo.Cartpole
using StaticArrays, LinearAlgebra
using RobotDynamics
const RD = RobotDynamics

# Use the Cartpole model from RobotZoo
model = Cartpole()
n,m = RD.dims(model)

# Define model discretization
N = 101
tf = 5.
dt = tf/(N-1)

# Define initial and final conditions
x0 = @SVector zeros(n)
xf = @SVector [0, pi, 0, 0]  # i.e. swing up

# Set up
Q = 1.0e-2*Diagonal(@SVector ones(n)) * dt
Qf = 100.0*Diagonal(@SVector ones(n))
R = 1.0e-1*Diagonal(@SVector ones(m)) * dt
obj = LQRObjective(Q,R,Qf,xf,N)

# Add constraints
conSet = ConstraintList(n,m,N)
u_bnd = 3.0
bnd = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd)
goal = GoalConstraint(xf)
add_constraint!(conSet, bnd, 1:N-1)
add_constraint!(conSet, goal, N)

# Initialization
u0 = @SVector fill(0.01,m)
U0 = [u0 for k = 1:N-1]

# Define problem
prob = Problem(model, obj, x0, tf, xf=xf, constraints=conSet)
initial_controls!(prob, U0)

# Solve with ALTRO
opts = SolverOptions(
    cost_tolerance_intermediate=1e-2,
    penalty_scaling=10.,
    penalty_initial=1.0
)
altro = ALTROSolver(prob, opts)
solve!(altro)

# Get some info on the solve
max_violation(altro)  # 5.896e-7
cost(altro)           # 1.539
iterations(altro)     # 44

# Extract the solution
X = states(altro)
U = controls(altro)

# Extract the solver statistics
stats = Altro.stats(altro)   # alternatively, solver.stats
stats.iterations             # 44, equivalent to iterations(solver)
stats.iterations_outer       # 4 (number of Augmented Lagrangian iterations)
stats.iterations_pn          # 1 (number of projected newton iterations)
stats.cost[end]              # terminal cost
stats.c_max[end]             # terminal constraint satisfaction
stats.gradient[end]          # terminal gradient of the Lagrangian
dstats = Dict(stats)         # get the per-iteration stats as a dictionary (can be converted to DataFrame)


text:
using DifferentialEquations, Test

function f(du,u,p,t)
  (x, v) = u
  du[1] = v
  du[2] = -x
end

function bc!(resid,sol,p,t)
  resid[1] = sol[1][1]
  resid[2] = sol[end][1] - 1
end

tspan = (0.,100.)
u0 = [0.,1.]
bvp = BVProblem(f, bc!, u0, tspan)
resid_f = Array{Float64}(undef, 2)
sol = solve(bvp, Shooting(Tsit5()))
sol2 = solve(bvp)

@test sol2.alg == Tsit5()
@test all(sol.u .== sol2.u)

max_stars_repo_name: ChrisRackauckas/DiffEq

max_issues_repo_name: ChrisRackauckas/DiffEq


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 66368ab4-ff5b-11ea-11f4-ad34524aa659
using Pkg, DrWatson

# ╔═╡ 6636c524-ff5b-11ea-1cdc-458faf68f551
begin
	@quickactivate "StatisticalRethinkingDHMC"
  using LogDensityProblems, DynamicHMC
  using Parameters, TransformVariables
	using StatisticalRethinking
end

# ╔═╡ 47f3e088-ff5b-11ea-033f-573e13956713
md"## Clip-09-03d.jl"

# ╔═╡ 66374e18-ff5b-11ea-3e4f-ed01da7ccf4f
md"### snippet 9.3"

# ╔═╡ 663dfcfe-ff5b-11ea-30f3-d7ce1a69f351
md"##### Construct the logdensity problem."

# ╔═╡ 664fa592-ff5b-11ea-1852-13135de62837
struct clip_9_3_model{TY <: AbstractVector, TX <: AbstractVector}
    "Observations."
    y::TY
    "Covariate"
    x::TX
end

# ╔═╡ 66506f7e-ff5b-11ea-2b80-0d40b56d40df
md"##### Make the type callable with the parameters *as a single argument*."

# ╔═╡ 665ba830-ff5b-11ea-3483-07fccf403fe4
function (problem:: clip_9_3_model)(θ)
    @unpack y, x, = problem    # extract the data
    @unpack muy, mux = θ     # works on the named tuple too
    ll = 0.0
    ll += loglikelihood(Normal(mux, 1), x)
    ll += loglikelihood(Normal(muy, 1), y)
    ll += logpdf(Normal(0, 1), mux)
    ll += logpdf(Normal(0, 1), muy)
    ll
end

# ╔═╡ 2c21e464-ff5a-11ea-1786-81348a9bd7c3
md"##### Instantiate the model with data and inits."

# ╔═╡ 814aae02-ff5b-11ea-2551-4ddbc8cc7290
begin
	Random.seed!(1234591)
	N = 100
	x = rand(Normal(0, 1), N)
	y = rand(Normal(0, 1), N)

	p = clip_9_3_model(y, x)
	θ = (muy = 0.0, mux=0.0)
	p(θ)
end

# ╔═╡ 814aefd4-ff5b-11ea-17b7-094af36740c3
md"##### Write a function to return properly dimensioned transformation."

# ╔═╡ 814b731e-ff5b-11ea-2fc4-9d3b3a48e82a
problem_transformation(p::clip_9_3_model) = as((muy = asℝ, mux = asℝ))

# ╔═╡ 815def08-ff5b-11ea-1464-fd9ed7306be7
md"##### Wrap the problem with a transformation, then use Flux for the gradient."

# ╔═╡ 815eade4-ff5b-11ea-3d15-bdaebd138d89
P = TransformedLogDensity(problem_transformation(p), p)

# ╔═╡ 816ac336-ff5b-11ea-0ba8-09a46f99b9b9
∇P = ADgradient(:ForwardDiff, P)

# ╔═╡ 816b7f74-ff5b-11ea-0287-f1c6e72dd350
md"##### Generate and show fig 9.3."

# ╔═╡ 8171f0fc-ff5b-11ea-3b97-2f07c9a9acdc
function generate_n_samples(model, grad;
  epsilon = 0.03,						# Step size
  L = 11,								# No of leapfrog steps
  pr = 0.3,								# Plot range
  n_samples = 4,						# No of samples
  q = [-0.1, 0.2])						# Initial position

  fig = plot(ylab="muy", xlab="mux", xlim=(-pr, pr), ylim=(-pr, pr), leg=false)

  for i in 1:n_samples
    q, ptraj, qtraj, accept, dH = HMC(model, grad, 0.03, 11, q)
    if n_samples < 10
      for j in 1:L
        K0 = sum(ptraj[j, :] .^ 2) / 2
        plot!(fig, [qtraj[j, 1], qtraj[j+1, 1]], [qtraj[j, 2], qtraj[j+1, 2]], leg=false)
        scatter!(fig, [(qtraj[j, 1], qtraj[j, 2])])
      end
    end
  end
  fig
end

# ╔═╡ 817f7d8c-ff5b-11ea-293d-2b2aad2823b2
generate_n_samples(p, ∇P; pr=0.5)

# ╔═╡ 818812e0-ff5b-11ea-361c-21dfc2292fd3
md"## End of clip-09-03d.jl"

# ╔═╡ Cell order:
# ╟─47f3e088-ff5b-11ea-033f-573e13956713
# ╠═66368ab4-ff5b-11ea-11f4-ad34524aa659
# ╠═6636c524-ff5b-11ea-1cdc-458faf68f551
# ╟─66374e18-ff5b-11ea-3e4f-ed01da7ccf4f
# ╟─663dfcfe-ff5b-11ea-30f3-d7ce1a69f351
# ╠═664fa592-ff5b-11ea-1852-13135de62837
# ╟─66506f7e-ff5b-11ea-2b80-0d40b56d40df
# ╠═665ba830-ff5b-11ea-3483-07fccf403fe4
# ╟─2c21e464-ff5a-11ea-1786-81348a9bd7c3
# ╠═814aae02-ff5b-11ea-2551-4ddbc8cc7290
# ╟─814aefd4-ff5b-11ea-17b7-094af36740c3
# ╠═814b731e-ff5b-11ea-2fc4-9d3b3a48e82a
# ╟─815def08-ff5b-11ea-1464-fd9ed7306be7
# ╠═815eade4-ff5b-11ea-3d15-bdaebd138d89
# ╠═816ac336-ff5b-11ea-0ba8-09a46f99b9b9
# ╟─816b7f74-ff5b-11ea-0287-f1c6e72dd350
# ╠═8171f0fc-ff5b-11ea-3b97-2f07c9a9acdc
# ╠═817f7d8c-ff5b-11ea-293d-2b2aad2823b2
# ╟─818812e0-ff5b-11ea-361c-21dfc2292fd3


text:
mutable struct AutoAbstolAffect{T}
  curmax::T
end
# Now make `affect!` for this:
function (p::AutoAbstolAffect)(integrator)
  if typeof(p.curmax) <: AbstractArray
    @. p.curmax = max(p.curmax,abs(integrator.u))
  else
    p.curmax = max(p.curmax,maximum(abs.(integrator.u)))
  end

  if typeof(integrator.opts.abstol) <: AbstractArray
    integrator.opts.abstol .= p.curmax .* integrator.opts.reltol
  else
    integrator.opts.abstol = p.curmax .* integrator.opts.reltol
  end

  u_modified!(integrator,false)
end

function AutoAbstol_initialize(cb,u,t,integrator)
  if cb.affect!.curmax == zero(integrator.opts.abstol)
    cb.affect!.curmax = integrator.opts.abstol
  end
  u_modified!(integrator,false)
end

function AutoAbstol(save=true;init_curmax=0.0)
  affect! = AutoAbstolAffect(abs.(init_curmax))
  condtion = (u,t,integrator) -> true
  save_positions = (save,false)
  DiscreteCallback(condtion,affect!;
                   initialize = AutoAbstol_initialize,
                   save_positions=save_positions)
end

export AutoAbstol


text:
push!(LOAD_PATH, "../src/")

using Documenter, CEDICT

makedocs(
    sitename="CEDICT.jl Documentation",
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    ),
    modules=[CEDICT],
    pages=[
        "Home" => "index.md",
        "API Reference" => [
            "Loading Dictionaries" => "api_dictionaries.md",
            "Searching in Dictionaries" => "api_searching.md",
            "Convenience Functions" => "api_convenience.md"
        ]
    ]
)

deploydocs(
    repo="github.com/tmthyln/CEDICT.jl.git",
    devbranch="main",
    devurl="latest"
    )


text:
@testset "Data model" begin
    @testset "Type dispatch" begin
        function thingwithstring(::GenserStringValue) true end

        @test thingwithstring(GenserString("OK"))
        intv = GenserInt8(1)
        @test_throws MethodError thingwithstring(intv)
    end

    @testset "Tag dispatch" begin
        function thingwithuri(::GenserDataType{Genser.uri}) true end

        @test thingwithuri(GenserURI("OK"))
        strv = GenserString("OK")
        @test_throws MethodError thingwithuri(strv)
    end

    @testset "Tag" begin
        @test Genser.tag(GenserURI) == Genser.uri
        @test Genser.tag(GenserNull()) == Genser.null
        @test Genser.tag(GenserNull) == Genser.null
        @test Genser.tag(GenserSequence{GenserNull}) == Genser.sequence
    end
end


text:
a = 1
for b = 2:20
  global a
  a = Int(a*b / gcd(a, b))
end
println(a)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
export Contour2D, points

"""
    Contour2D

Represents a contour of a two-dimensional amoeba.
"""
struct Contour2D
    ps::Vector{SVector{2, Float64}}
    grid::Grid2D
end

"""
    points(C::Contour2D)

The sample points describing `C`.
"""
points(C::Contour2D) = C.ps

function Base.show(io::IO, A::Contour2D)
    println(io, "Contour2D:")
    println(io, " * domain: ", limits(A.grid))
    println(io, " * number of points: ", length(points(A)))
end

function Base.show(io::IO, mime::MIME"text/html", A::Contour2D)
    if inIJulia() && plotsdefined()
        show(io, mime, Main.Plots.plot(A))
    else
        show(io, A)
    end
end

@recipe function plot(C::Contour2D)
    xmin, xmax, ymin, ymax = limits(C.grid)
    @series begin

        if !haskey(plotattributes, :xlims)
            xlims --> (xmin, xmax)
        end
        if !haskey(plotattributes, :ylims)
            ylims --> (ymin, ymax)
        end
        aspect_ratio --> :equal
        # grid --> false
        legend := false
        st := :scatter
        #markershape := :rect
        markercolor --> Colors.colorant"#884EA0"
        markerstrokewidth := 0.0
        markersize --> 2.0

        xs, ys = unzip(C.ps)
        xs, ys
    end

end

function _contour2D(F::ContourFiber2D, grid::Grid2D,
    samples_off_axis,
    options::MembershipTestOptions)

    ntries = options.ntries
    iterations = options.maxiters
    tol = options.tol

    fix_axis!(F, :w1)
    Fs = [deepcopy(F) for _=2:Threads.nthreads()]
    push!(Fs, F)

    threads_ps = [Vector{SVector{2, Float64}}() for _=1:Threads.nthreads()]
    threads_ws = [Vector{Float64}() for _=1:Threads.nthreads()]

    xmin, xmax, ymin, ymax = limits(grid)
    # first pass
    w1s = xcoordinates(grid)
    w2s = range(ymin, stop=ymax, length=samples_off_axis)

    force_exit_y(x) = x[3] > ymax || x[3] < ymin

    Threads.@threads for k=1:length(w1s)
        G = Fs[Threads.threadid()]
        w1_ys = threads_ws[Threads.threadid()]
        ps = threads_ps[Threads.threadid()]
        w1 = w1s[k]
        update_fiber!(G, w1)

        for w2 in w2s
            k = 0
            while k < ntries
                k += 1
                x = SVector(2π*rand(), 2π*rand(), w2)
                res = newton(G, x, iterations, tol, force_exit_y)
                if res.converged
                    y = res.solution[3]
                    p = SVector(w1, y)
                    if isunique(w1_ys, y, step(w1s) * 0.5)
                        push!(ps, p)
                        push!(w1_ys, y)
                    end
                end
            end
        end
        empty!(w1_ys)
    end

    # second pass
    for F in Fs
        fix_axis!(F, :w2)
    end
    w1s = range(xcoordinates(grid)[1], stop=xcoordinates(grid)[end], length=samples_off_axis)
    w2s = ycoordinates(grid)
    # w2_xs = Vector{Float64}()
    force_exit_x(x) = x[3] > xmax || x[3] < xmin

    Threads.@threads for k=1:length(w2s)
        G = Fs[Threads.threadid()]
        w2_xs = threads_ws[Threads.threadid()]
        ps = threads_ps[Threads.threadid()]
        w2 = w2s[k]
        update_fiber!(G, w2)
        for w1 in w1s
            k = 0
            while k < ntries
                k += 1
                x = SVector(2π*rand(), 2π*rand(), w1)
                res = newton(G, x, iterations, tol, force_exit_x)
                if res.converged
                    x = res.solution[3]
                    p = SVector(x, w2)
                    if isunique(w2_xs, x, step(w2s) * 0.5)
                        push!(ps, p)
                        push!(w2_xs, x)
                    end
                end
            end
        end
        empty!(w2_xs)
    end

    ps = Vector{SVector{2, Float64}}()
    for i=1:Threads.nthreads()
        append!(ps, threads_ps[i])
    end
    tol = 0.5 * min(step(xcoordinates(grid)), step(ycoordinates(grid)))
    ps = unique_points(ps, tol)

    Contour2D(ps, grid)
end

function _contour2D_grid(F::ContourFiber2D, grid::Grid2D,
    every_main_axis, samples_off_axis, options::MembershipTestOptions)

    ntries = options.ntries
    iterations = options.maxiters
    tol = options.tol

    fix_axis!(F, :w1)
    Fs = [deepcopy(F) for _=2:Threads.nthreads()]
    pushfirst!(Fs, F)

    B = Bitmap2D(grid)

    xmin, xmax, ymin, ymax = limits(grid)

    # FIRST PASS
    w1s = xcoordinates(grid)
    w2s = central_range(ymin, ymax, samples_off_axis)

    force_exit_y(x) = x[3] > ymax || x[3] < ymin

    Threads.@threads for k=1:every_main_axis:length(w1s)
        G = Fs[Threads.threadid()]
        w1 = w1s[k]
        update_fiber!(G, w1)
        px = xgridpoint(grid, w1)
        for w2 in w2s
            k = 0
            while k < ntries
                k += 1
                x = SVector(2π*rand(), 2π*rand(), w2)
                res = newton(G, x, iterations, tol, force_exit_y)
                if res.converged
                    py = ygridpoint(grid, res.solution[3])
                    safe_setindex!(B, true, py, px)
                end
            end
        end
    end

    # SECOND PASS
    for F in Fs
        fix_axis!(F, :w2)
    end
    w1s = central_range(xmin, xmax, samples_off_axis)
    w2s = ycoordinates(grid)

    force_exit_x(x) = x[3] > xmax || x[3] < ymin

    Threads.@threads for k=1:every_main_axis:length(w2s)
        G = Fs[Threads.threadid()]
        w2 = w2s[k]
        update_fiber!(G, w2)
        py = ygridpoint(grid, w2)
        for w1 in w1s
            k = 0
            while k < ntries
                k += 1
                x = SVector(2π*rand(), 2π*rand(), w1)
                res = newton(G, x, iterations, tol, force_exit_x)
                if res.converged
                    px = xgridpoint(grid, res.solution[3])
                    safe_setindex!(B, true, py, px)
                end
            end
        end
    end

    B
end

function isunique(vec::Vector{Float64}, x::Float64, tol::Float64)
    for y in vec
        if abs(x - y) < tol
            return false
        end
    end
    return true
end

# function unique_tol(vec::AbstractVector, tol::Float64)
#     out =
#     for y in vec
#         if abs(x - y) < tol
#             return false
#         end
#     end
#     return true
# end


text:
function evalwithx(expr::Expr, a, b)
    a = eval(quote let x = $a; return $expr end end)
    b = eval(quote let x = $b; return $expr end end)
    return b - a
end
evalwithx(expr::AbstractString, a, b) = evalwithx(parse(expr), a, b)

evalwithx(:(2 ^ x), 3, 5)
evalwithx("2 ^ x", 3, 5)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function whiteNoiseGenerator(T, interval, sigma, seed=123)
	Random.seed!(seed)
	d = Normal(0,sigma^2/interval)
	return rand(d, Int(T/interval))
end


text:
# DOCOPT faster than argparse, but a little sensitive to spacing
doc = """We use this to find an MPPI trajectory. UW CSE Motion Control Lab.

Usage:
   PISQTrain.jl [(<experiment> <name>)]
   PISQTrain.jl -h | --help

Options:
   -h --help   Show this screen.

"""

using DocOpt
opt = docopt(doc, ARGS)

@everywhere push!(LOAD_PATH, "./")

using Common: plotlines

prefix="/tmp/"
modeldir="$(pwd())/models"
meshfile = "$(modeldir)/meshes/" # hack for now...

expfile=opt["<experiment>"]
expname=opt["<name>"]
myseed = 12345

@eval @everywhere const SEED = $myseed
@everywhere push!(LOAD_PATH, "$(pwd())")

println("Experiment called $(expname)")
println("Loading experiment parameters from $(expfile)")

################################################### start work by loading files
using mjWrap
using MPPI
using ExpFunctions
using LearningStrategies

@eval @everywhere const mjw = mjWrap # scoping issues...

const METHOD = :MPPI
include("$(pwd())/$(expfile)") # LOADING TASK FILE with PARAMETERS
@assert isdefined(:mppi)
@assert isdefined(:mppi_specs)
@assert isdefined(:myfuncs)

# plot stoc & mean every other iter
evaluatepolicy = IterFunction((model,i)->push!(model.trace[:evalscore], myfuncs.evaluate(model.mjsys, model.samples)))

const NITER = 100
plotR = IterFunction(NITER, (model,i)->plotlines(i,"Costs",
                                                 (model.trace[:meanR],"Min"),
                                                 (model.trace[:stocR],"Avg")))

plotEval = IterFunction(NITER, (model,i)->plotlines(i,"Evaluation",
                                                    (model.trace[:evalscore],"")))

plotCTRL = IterFunction(NITER, (model,i)->plotlines(i,"Ctrl",
                                                    model.samples.ctrl[:,:,1]))

plotState = IterFunction(NITER, (model,i)->plotlines(i,"State",
                                                     model.traj.state[1:3,1:i,1])) # X axis plot for now

# meta strategies: order matters!!
mppi_strat = strategy(ControlExp(prefix, expfile,
                                 expname, meshfile), # setup some things, hook for save
                      MPPIRollout(myfuncs),
                      mppi_specs,                  # main loop

                      #evaluatepolicy,             # hooks

                      plotR,

                      #plotEval,
                      plotCTRL,

                      #plotState,
                      MaxIter(mppi.niter)) #Verbose(MaxIter(mppi.niter)))

@time learn!(mppi, mppi_strat)

########################## done


text:
"Convert the indices in a graph to 32-bit ints.  This takes less storage, but does not speed up much"
shortIntGraph(a::SparseMatrixCSC) = SparseMatrixCSC{Float64,Int32}(convert(Int32,a.m), convert(Int32,a.n), convert(Array{Int32,1},a.colptr), convert(Array{Int32,1},a.rowval), a.nzval)

"Convert the nonzero entries in a graph to Float64"
floatGraph(a::SparseMatrixCSC) = SparseMatrixCSC{Float64,Int64}(a.m, a.n, a.colptr, a.rowval, convert(Array{Float64,1},a.nzval))

"""Create a Laplacian matrix from an adjacency matrix.
We might want to do this differently, say by enforcing symmetry"""
lap(a) = spdiagm(a*ones(size(a)[1])) - a

"""Create an adjacency matrix and a diagonal vector from a Laplacian with added diagonal weights"""
function adj{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti})
  a = -(2 * la - triu(la) - tril(la))
  d = diag(la) - diag(lap(a))
  return a,d
end

"""Create a new graph in that is the same as the original, but with all edge weights 1"""
function unweight{Tval,Tind}(ain::SparseMatrixCSC{Tval,Tind})
    a = copy(ain)
    m = length(a.nzval)
    for i in 1:m
        a.nzval[i] = 1
    end
    return a
end # unweight

"""Change the weight of every edge in a to 1"""
function unweight!{Tval,Tind}(a::SparseMatrixCSC{Tval,Tind})
    m = length(a.nzval)
    for i in 1:m
        a.nzval[i] = 1
    end
end # unweight

"""Create a new graph that is the same as the original, but with f applied to each nonzero entry of a. For example, to make the weight of every edge uniform in [0,1], we could write

~~~julia
b = mapweight(a, x->rand(1)[1])
~~~
"""
function mapweight{Tval,Tind}(a::SparseMatrixCSC{Tval,Tind},f)
    b = triu(a,1)
    b.nzval = map(f, b.nzval)
    b = b + b'

    return b

end # mapweight

"""Put a uniform [0,1] weight on every edge.  This is an example of how to use mapweight."""
uniformWeight{Tval,Tind}(a::SparseMatrixCSC{Tval,Tind}) = mapweight(a,x->rand(1)[1])

"""Set the weight of every edge to 1"""
function uniformWeight!(mat::SparseMatrixCSC)
    mat.nzval = ones(length(mat.nzval))
end

"""The Cartesian product of two graphs.  When applied to two paths, it gives a grid."""
function productGraph(a0::SparseMatrixCSC, a1::SparseMatrixCSC)
  n0 = size(a0)[1]
  n1 = size(a1)[1]
  a = kron(speye(n0),a1) + kron(a0,speye(n1));

end # productGraph

"""The signed edge-vertex adjacency matrix"""
function edgeVertexMat(mat::SparseMatrixCSC)
    (ai,aj) = findnz(triu(mat,1))
    m = length(ai)
    n = size(mat)[1]
    return sparse(collect(1:m),ai,1.0,m,n) - sparse(collect(1:m),aj,1.0,m,n)
end

"""Create a new graph from the old, but keeping edge edge with probability `p`"""
function subsampleEdges(a::SparseMatrixCSC{Float64,Int64}, p::Float64)
  (ai, aj, av) = findnz(triu(a))
  n = size(a)[1]
  mask = map(x -> convert(Float64,x < p), rand(length(ai)))
  as = sparse(ai,aj,mask.*av,n,n)
  as = as + as'
  return as

end # subsampleEdges

"""Creats a 2-lift of a.  `flip` is a boolean indicating which edges cross"""
function twoLift(a, flip::AbstractArray{Bool,1})
    (ai,aj,av) = findnz(triu(a))
    m = length(ai)
    #flip = rand(false:true,m)
    n = size(a)[1]
    a0 = sparse(ai[flip],aj[flip],1,n,n)
    a1 = sparse(ai[!flip],aj[!flip],1,n,n)
    a00 = a0 + a0'
    a11 = a1 + a1'
    return [a00 a11; a11 a00]
end

twoLift(a) = twoLift(a,rand(false:true,div(nnz(a),2)))

twoLift(a, k::Integer) = twoLift(a,randperm(div(nnz(a),2)) .> k)

"""
 Create a disjoint union of graphs a and b,
 and then put k random edges between them
"""
function joinGraphs{Tval,Tind}(a::SparseMatrixCSC{Tval,Tind}, b::SparseMatrixCSC{Tval,Tind}, k::Integer)
    na = size(a)[1]
    nb = size(b)[1]

    (ai,aj,av) = findnz(a)
    (bi,bj,bv) = findnz(b)
    bi = bi + na
    bj = bj + na

    ji = rand(1:na,k)
    jj = rand(1:nb,k)+na

    ab = sparse([ai;bi;ji;jj],[aj;bj;jj;ji],[av;bv;ones(Tval,2*k)],na+nb,na+nb)
end

"""
 Create a disjoint union of graphs a and b,
  with no edges between them.
"""
disjoin(a,b) = joinGraphs(a,b,0)

"""Plots graph gr with coordinates (x,y)"""
function plotGraph(gr,x,y,color=[0,0,1];dots=true,setaxis=true,number=false)
  (ai,aj,av) = findnz(triu(gr))
  arx = [x[ai]';x[aj]';NaN*ones(length(ai))']
  ary = [y[ai]';y[aj]';NaN*ones(length(ai))']
  p = plot(arx[:],ary[:],color=color)

  if dots
    plot(x,y,color=color,marker="o",linestyle="none")
  end #if

  if number
    for i in 1:length(x)
      annotate(i, xy = [x[i]; y[i]])
    end
  end

  axis("off")

  if setaxis
  ax = axes()

  minx = minimum(x)
  maxx = maximum(x)
  miny = minimum(y)
  maxy = maximum(y)
  delx = maxx - minx
  dely = maxy - miny
  ax[:set_ylim]([miny - dely/20, maxy + dely/20])
  ax[:set_xlim]([minx - delx/20, maxx + delx/20])
  end

  return p
end # plotGraph

"""Computes spectral coordinates, and then uses plotGraph to draw"""
function spectralDrawing(a)

    x, y = spectralCoords(a)
    plotGraph(a,x,y)

end # spectralDrawing

"""Computes the spectral coordinates of a graph"""
function spectralCoords(a)

    E = eigs(lap(a), nev = 3, which=:SR)
    V = E[2]
    return V[:,2], V[:,3]

end # spectralCoords

"""Creates a unit vector of length n from a given set of integers, with weights based on the number of occurences"""
function toUnitVector(a::Array{Int64,1}, n)

  v = zeros(n)

  for i in a
    v[i] = v[i] + 1
  end

  sum = 0
  for i in 1:length(v)
    sum = sum + v[i] * v[i]
  end

  for i in 1:length(v)
    v[i] = v[i] / sqrt(sum)
  end

  return v

end # toUnitVector

"Returns the diagonal matrix(as a sparse matrix) of a graph"
function diagmat{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti})

  dw = zeros(Tv, G.n)

  u,v,w = findnz(G)

  for i in 1:length(u)
    dw[v[i]] = dw[v[i]] + w[i]
  end

  return sparse(collect(1:G.n), collect(1:G.n), dw)

end # diagmat

"""
 Constructs a generalized necklace graph starting with two graphs A and H. The
resulting new graph will be constructed by expanding each vertex in H to an
instance of A. k random edges will be generated between components. Thus, the
resulting graph may have weighted edges.
"""
function generalizedNecklace{Tv, Ti}(A::SparseMatrixCSC{Tv, Ti}, H::SparseMatrixCSC, k::Int64)
  a = findnz(A)
  h = findnz(H)

  # these are square matrices
  n = A.n
  m = H.n

  newI = Ti[]
  newJ = Ti[]
  newW = Tv[]

  # duplicate the vertices in A so that each vertex in H corresponds to a copy of A
  for i in 1:m
    newI = append!(newI, a[1] + n * (i - 1))
    newJ = append!(newJ, a[2] + n * (i - 1))
    newW = append!(newW, a[3])
  end

  # for each edge in H, add k random edges between two corresponding components
  # multiedges will be concatenated to a single edge with higher cost
  for i in 1:length(h[1])
    u = h[1][i]
    v = h[2][i]

    if (u < v)
      #component x is from 1 + (x - 1) * n to n + (x - 1) * n
      for edgeToAdd in 1:k
        newU = rand(1:n) + n * (u - 1)
        newV = rand(1:n) + n * (v - 1)
        append!(newI, [newU, newV])
        append!(newJ, [newV, newU])
        append!(newW, [1, 1])
      end
    end
  end

  return sparse(newI, newJ, newW)
end # generalizedNecklace


text:
module ContractionSequenceOptimization

export optimal_contraction_sequence, contraction_cost

include("utils.jl")
include("three_tensors.jl")
include("breadth_first_constructive.jl")
include("contraction_cost.jl")

#
# The integer type of the dimensions and costs.
# Needs to be large to avoid overflow, but larger
# integer types are a little bit slower.
#
# For very large dimensions, one could try:
#
# https://github.com/rfourquet/BitIntegers.jl
#
# however in that limit it would be best to use
# symbolic dimensions.
#
const DimT = UInt128

"""
    optimal_contraction_sequence(T)

Returns a contraction sequence for contracting the tensors `T`. The sequence is generally optimal (currently, outer product contractions are skipped, but some optimal sequences require outer product contractions).
"""
function optimal_contraction_sequence(T)
  if length(T) == 1
    return Any[1]
  elseif length(T) == 2
    return Any[1, 2]
  elseif length(T) == 3
    return optimal_contraction_sequence(T[1], T[2], T[3])
  end
  return breadth_first_constructive(T)
end

end # module ContractionSequenceOptimization