text:
# Written by Jack Wisdom in June and July, 2015.
# David M. Hernandez helped test and improve it.
# Translated to Julia by Eric Agol August 2018.

#include "universal.h"

int solve_universal_newton(double kc, double r0, double beta, double b, double eta, double zeta, double h,
			   double *X, double *S2, double *C2)
{
  double x;
  double g1, g2, g3, arg;
  double s2, c2, xnew;
  int count;
  double cc;
  double adx;

  xnew = *X;

  count = 0;
  do {
    x = xnew;
    arg = b*x/2.0;
    s2 = sin(arg);
    c2 = cos(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/beta;
    g3 = (x - g1)/beta;
    cc = eta*g1 + zeta*g2;
    xnew = (h + (x*cc - (eta*g2 + zeta*g3)))/(r0 + cc);
    if(count++ > 10) {
      return(FAILURE);
    }
  } while(fabs((x-xnew)/xnew) > 1.e-8);

  /* compute the output */
  x = xnew;
  arg = b*x/2.0;
  s2 = sin(arg);
  c2 = cos(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

int solve_universal_laguerre(double kc, double r0, double beta, double b, double eta, double zeta, double h,
			     double *X, double *S2, double *C2)
{
  double x;
  double g1, g2, g3;
  double arg, s2, c2, xnew, dx;
  double f, fp, fpp, g0;
  int count;
  double cc;

  xnew = *X;

  count = 0;
  do {
    double c5=5.0, c16 = 16.0, c20 = 20.0;

    x = xnew;
    arg = b*x/2.0;
    s2 = sin(arg);
    c2 = cos(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/beta;
    g3 = (x - g1)/beta;
    f = r0*x + eta*g2 + zeta*g3 - h;
    fp = r0 + eta*g1 + zeta*g2;
    g0 = 1.0 - beta*g2;
    fpp = eta*g0 + zeta*g1;
    dx = -c5*f/(fp + sqrt(fabs(c16*fp*fp - c20*f*fpp)));
    xnew = x + dx;
    if(count++ > 10) {
      return(FAILURE);
    }
  } while(fabs(dx) > 2.e-7*fabs(xnew));

  /* refine with a last newton */
  x = xnew;
  arg = b*x/2.0;
  s2 = sin(arg);
  c2 = cos(arg);
  g1 = 2.0*s2*c2/b;
  g2 = 2.0*s2*s2/beta;
  g3 = (x - g1)/beta;
  cc = eta*g1 + zeta*g2;
  xnew = (h + (x*cc - (eta*g2 + zeta*g3)))/(r0 + cc);

  x = xnew;
  arg = b*x/2.0;
  s2 = sin(arg);
  c2 = cos(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

int solve_universal_bisection(double kc, double r0, double beta, double b, double eta, double zeta, double h,
				     double *X, double *S2, double *C2)
{
  double x;
  double g1, g2, g3;
  double arg, s2, c2;
  double f;
  double X_min, X_max, X_per_period, invperiod;
  int count;
  double xnew;
  double err;

  xnew = *X;

  err = 1.e-9*fabs(xnew);

  /* bisection limits due to Rein et al. */
  invperiod = b*beta/(2.*M_PI*kc);
  X_per_period = 2.*M_PI/b;
  X_min = X_per_period * floor(h*invperiod);
  X_max = X_min + X_per_period;
  xnew = (X_max + X_min)/2.;

  count = 0;
  do {
    x = xnew;
    arg = b*x/2.0;
    s2 = sin(arg);
    c2 = cos(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/beta;
    g3 = (x - g1)/beta;
    f = r0*x + eta*g2 + zeta*g3 - h;
    if (f>=0.){
      X_max = x;
    }else{
      X_min = x;
    }
    xnew = (X_max + X_min)/2.;
    if(count++ > 100) return(FAILURE);

  } while(fabs(x - xnew) > err);

  x = xnew;

  arg = b*x/2.0;
  s2 = sin(arg);
  c2 = cos(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

double sign(double x)
{
  if(x > 0.0) {
    return(1.0);
  } else if(x < 0.0) {
    return(-1.0);
  } else {
    return(0.0);
  }
}

function cubic1(a::T, b::T, c::T) where{T <: Real}
{
  double Q, R;
  double theta, A, B, x1;

  Q = (a*a - 3.0*b)/9.0;
  R = (2.0*a*a*a - 9.0*a*b + 27.0*c)/54.0;
  if(R*R < Q*Q*Q)
    theta = acos(R/sqrt(Q*Q*Q));

    println("three cubic roots %.16le %.16le %.16le\n",
	   -2.0*sqrt(Q)*cos(theta/3.0) - a/3.0,
	   -2.0*sqrt(Q)*cos((theta + 2.0*M_PI)/3.0) - a/3.0,
	   -2.0*sqrt(Q)*cos((theta - 2.0*M_PI)/3.0) - a/3.0);
    exit(-1)

    return x1
  else
    A = -sign(R)*pow(fabs(R) + sqrt(R*R - Q*Q*Q), 1./3.);
    if(A == 0.0) {
      B = 0.0;
    } else {
      B = Q/A;
    }
    x1 = (A + B) - a/3.0;

    return(x1);
  end
end

int solve_universal_parabolic(double kc, double r0, double beta, double b, double eta, double zeta, double h,
			      double *X, double *S2, double *C2)
{
  double x;
  double s2, c2;

  b = 0.0;

  s2 = 0.0;
  c2 = 1.0;
  /*
    g1 = x;
    g2 = x*x/2.0;
    g3 = x*x*x/6.0;
    g = eta*g1 + zeta*g2;
   */
  /*   f = r0*x + eta*g2 + zeta*g3 - h; */
  /*   this is just a cubic equation */
  x = cubic1(3.0*eta/zeta, 6.0*r0/zeta, -6.0*h/zeta);

  s2 = 0.0;
  c2 = 1.0;
  /* g1 = 2.0*s2*c2/b; */
  /* g2 = 2.0*s2*s2/beta; */

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

int solve_universal_hyperbolic_newton(double kc, double r0, double minus_beta, double b, double eta, double zeta, double h,
				      double *X, double *S2, double *C2)
{
  double x;
  double g1, g2, g3;
  double arg, s2, c2, xnew;
  int count;
  double g;

  xnew = *X;

  count = 0;
  do {
    x = xnew;
    arg = b*x/2.0;
    if(fabs(arg)>200.0) return(FAILURE);
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    g = eta*g1 + zeta*g2;
    xnew = (x*g - eta*g2 - zeta*g3 + h)/(r0 + g);

    if(count++ > 10) return(FAILURE);

  } while(fabs(x - xnew) > 1.e-9*fabs(xnew));

  x = xnew;
  arg = b*x/2.0;
  s2 = sinh(arg);
  c2 = cosh(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

int solve_universal_hyperbolic_laguerre(double kc, double r0, double minus_beta, double b, double eta, double zeta, double h,
					double *X, double *S2, double *C2)
{
  double x;
  double g0, g1, g2, g3;
  double arg, s2, c2, xnew;
  int count;
  double f;

  xnew = *X;

  count = 0;
  do {
    double c5=5.0, c16 = 16.0, c20 = 20.0;
    double fp, fpp, dx;
    double den;

    x = xnew;
    arg = b*x/2.0;
    if(fabs(arg)>50.0) return(FAILURE);
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    f = r0*x + eta*g2 + zeta*g3 - h;
    fp = r0 + eta*g1 + zeta*g2;
    g0 = 1.0 + minus_beta*g2;
    fpp = eta*g0 + zeta*g1;
    den = (fp + sqrt(fabs(c16*fp*fp - c20*f*fpp)));
    if(den == 0.0) return(FAILURE);
    dx = -c5*f/den;
    xnew = x + dx;
    if(count++ > 20) return(FAILURE);

  } while(fabs(x - xnew) > 1.e-9*fabs(xnew));

  /* refine with a step of Newton */
  {
    double g;

    x = xnew;
    arg = b*x/2.0;
    if(fabs(arg)>200.0) return(FAILURE);
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    g = eta*g1 + zeta*g2;
    xnew = (x*g - eta*g2 - zeta*g3 + h)/(r0 + g);
  }

  x = xnew;
  arg = b*x/2.0;
  s2 = sinh(arg);
  c2 = cosh(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

int solve_universal_hyperbolic_bisection(double kc, double r0, double minus_beta, double b, double eta, double zeta, double h,
					 double *X, double *S2, double *C2)
{
  double x;
  double g1, g2, g3;
  double arg, s2, c2, xnew;
  int count;
  double X_min, X_max;
  double f;
  double err;

  xnew = *X;

  err = 1.e-10*fabs(xnew);

  X_min = 0.5*xnew;
  X_max = 10.0*xnew;

  {
    double fmin, fmax;
    x = X_min;
    arg = b*x/2.0;
    if(fabs(arg)>200.0) return(FAILURE);
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    fmin = r0*x + eta*g2 + zeta*g3 - h;

    x = X_max;
    arg = b*x/2.0;
    if(fabs(arg)>200.0) {
      x = 200.0/(b/2.0);
      arg = 200.0;
    }
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    fmax = r0*x + eta*g2 + zeta*g3 - h;

    if(fmin*fmax > 0.0) {
      return(FAILURE);
    }
  }

  count = 0;
  do {
    x = xnew;
    arg = b*x/2.0;
    if(fabs(arg)>200.0) return(FAILURE);
    s2 = sinh(arg);
    c2 = cosh(arg);
    g1 = 2.0*s2*c2/b;
    g2 = 2.0*s2*s2/minus_beta;
    g3 = -(x - g1)/minus_beta;
    f = r0*x + eta*g2 + zeta*g3 - h;
    if (f>=0.){
      X_max = x;
    }else{
      X_min = x;
    }
    xnew = (X_max + X_min)/2.;
    if(count++ > 100) return(FAILURE);

  } while(fabs(x - xnew) > err);

  x = xnew;
  arg = b*x/2.0;
  s2 = sinh(arg);
  c2 = cosh(arg);

  *X = x;
  *S2 = s2;
  *C2 = c2;

  return(SUCCESS);
}

void kepler_step(double kc, double dt, State *s0, State *s)
{
  void kepler_step_depth();
  double r0, v2, eta, beta, zeta, b;

  r0 = sqrt(s0->x*s0->x + s0->y*s0->y + s0->z*s0->z);
  v2 = s0->xd*s0->xd + s0->yd*s0->yd + s0->zd*s0->zd;
  eta = s0->x*s0->xd + s0->y*s0->yd + s0->z*s0->zd;
  beta = 2.0*kc/r0 - v2;
  zeta = kc - beta*r0;
  b = sqrt(fabs(beta));

  kepler_step_depth(kc, dt, beta, b, s0, s, 0, r0, v2, eta, zeta);
}

void kepler_step_depth(double kc, double dt, double beta, double b,
		       State *s0, State *s, int depth,
		       double r0, double v2, double eta, double zeta)
{
  int flag;
  State ss;
  int kepler_step_internal();

  if(depth > 30) {
    printf("kepler depth exceeded\n");
    exit(-1);
  }

  flag = kepler_step_internal(kc, dt, beta, b, s0, s, r0, v2, eta, zeta);

  if(flag == FAILURE) {
    kepler_step_depth(kc, dt/4.0, beta, b, s0, &ss, depth+1, r0, v2, eta, zeta);

    r0 = sqrt(ss.x*ss.x + ss.y*ss.y + ss.z*ss.z);
    v2 = ss.xd*ss.xd + ss.yd*ss.yd + ss.zd*ss.zd;
    eta = ss.x*ss.xd + ss.y*ss.yd + ss.z*ss.zd;
    zeta = kc - beta*r0;

    kepler_step_depth(kc, dt/4.0, beta, b, &ss, s, depth+1, r0, v2, eta, zeta);

    r0 = sqrt(s->x*s->x + s->y*s->y + s->z*s->z);
    v2 = s->xd*s->xd + s->yd*s->yd + s->zd*s->zd;
    eta = s->x*s->xd + s->y*s->yd + s->z*s->zd;
    zeta = kc - beta*r0;

    kepler_step_depth(kc, dt/4.0, beta, b, s, &ss, depth+1, r0, v2, eta, zeta);

    r0 = sqrt(ss.x*ss.x + ss.y*ss.y + ss.z*ss.z);
    v2 = ss.xd*ss.xd + ss.yd*ss.yd + ss.zd*ss.zd;
    eta = ss.x*ss.xd + ss.y*ss.yd + ss.z*ss.zd;
    zeta = kc - beta*r0;

    kepler_step_depth(kc, dt/4.0, beta, b, &ss, s, depth+1, r0, v2, eta, zeta);
  }
}

function new_guess(r0::T, eta::T, zeta::T, dt::T) where {T <: Real}
  nil = zero(T)
  if (zeta != zro)
    s = cubic1(3eta/zeta, 6r0/zeta, -6dt/zeta)
  elseif (eta != nil)
    reta = r0/eta
    disc = reta*reta + 2dt/eta
    if (disc >= nil)
      s = -reta + sqrt(disc)
    else
      s = dt/r0
    end
  else
    s = dt/r0
  end
  return s
end

int kepler_step_internal(double kc, double dt, double beta, double b,
			 State *s0, State *s,
			 double r0, double v2, double eta, double zeta)
{
  double r;
  double c2, s2;
  double G1, G2;
  double a, x;
  int flag;
  double c, ca, bsa;

  if(beta < 0.0) {
    double x0;

    x0 = new_guess(r0, eta, zeta, dt);

    x = x0;
    flag = solve_universal_hyperbolic_newton(kc, r0, -beta, b, eta, zeta, dt, &x, &s2, &c2);

    if(flag == FAILURE) {
      x = x0;
      flag = solve_universal_hyperbolic_laguerre(kc, r0, -beta, b, eta, zeta, dt, &x, &s2, &c2);
    }

    /*
    if(flag == FAILURE) {
      x = x0;
      flag = solve_universal_hyperbolic_bisection(kc, r0, -beta, b, eta, zeta, dt, &x, &s2, &c2);
    }
      */
    if(flag == FAILURE) {
      return(FAILURE);
    }

    a = kc/(-beta);
    G1 = 2.0*s2*c2/b;
    c = 2.0*s2*s2;
    G2 = c/(-beta);
    ca = c*a;
    r = r0 + eta*G1 + zeta*G2;
    bsa = (a/r)*(b/r0)*2.0*s2*c2;

  } else if(beta > 0.0) {
    double x0;
    double ff, fp;

    /* x0 = dt/r0; */

    x0 = dt/r0;
    ff = zeta*x0*x0*x0 + 3.0*eta*x0*x0;
    fp = 3.0*zeta*x0*x0 + 6.0*eta*x0 + 6.0*r0;
    x0 -= ff/fp;

    x = x0;
    flag = solve_universal_newton(kc, r0, beta, b, eta, zeta, dt, &x, &s2, &c2);
    if(flag == FAILURE) {
      x = x0;
      flag = solve_universal_laguerre(kc, r0, beta, b, eta, zeta, dt, &x, &s2, &c2);
    }
    /*
    if(flag == FAILURE) {
      x = x0;
      flag = solve_universal_bisection(kc, r0, beta, b, eta, zeta, dt, &x, &s2, &c2);
    }
      */
    if(flag == FAILURE) {
      return(FAILURE);
    }

    a = kc/beta;
    G1 = 2.0*s2*c2/b;
    c = 2.0*s2*s2;
    G2 = c/beta;
    ca = c*a;
    r = r0 + eta*G1 + zeta*G2;
    bsa = (a/r)*(b/r0)*2.0*s2*c2;

  } else {

    x = dt/r0;

    flag = solve_universal_parabolic(kc, r0, beta, b, eta, zeta, dt, &x, &s2, &c2);
    if(flag == FAILURE) {
      exit(-1);
    }

    G1 = x;
    G2 = x*x/2.0;
    ca = kc*G2;
    r = r0 + eta*G1 + zeta*G2;
    bsa = kc*x/(r*r0);

  }

  {
    double g, fdot, fhat, gdothat;

    /* f = 1.0 - (ca/r0); */
    fhat = -(ca/r0);
    g = eta*G2 + r0*G1;
    fdot = -bsa;
    /* gdot = 1.0 - (ca/r); */
    gdothat = -(ca/r);

    s->x = s0->x + (fhat*s0->x + g*s0->xd);
    s->y = s0->y + (fhat*s0->y + g*s0->yd);
    s->z = s0->z + (fhat*s0->z + g*s0->zd);
    s->xd = s0->xd + (fdot*s0->x + gdothat*s0->xd);
    s->yd = s0->yd + (fdot*s0->y + gdothat*s0->yd);
    s->zd = s0->zd + (fdot*s0->z + gdothat*s0->zd);
  }

  return(SUCCESS);
}


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofRepFAM
using CytofRepFAM, Random, BSON, Test, Distributions
=#

function init_state_const_data(; N=[300, 200, 100], J=8, K=4,
                               L=Dict(0=>5, 1=>3), KMCMC=K*2,
                               sig2=fill(0.1, length(N)),
                               mus=Dict(0=>collect(range(-5, length=L[0], stop=-1)),
                                        1=>collect(range(1, length=L[1], stop=5))),
                               LMCMC=Dict(0=>5, 1=>5), seed=-1,
                               allow_repeated_Z_columns=true)

  if seed >= 0
    Random.seed!(seed)
  end

  I = length(N)
  simdat = CytofRepFAM.Model.genData(J, N, K, L,
                                     mus=mus, sig2=sig2,
                                     sortLambda=true)
  d = CytofRepFAM.Model.Data(simdat[:y])
  simz = CytofRepFAM.Model.sim_fn_exp_decay_generator(1.0)
  deltaz_prior = TruncatedNormal(1.0, 0.1, 0.0, Inf)
  c = CytofRepFAM.Model.defaultConstants(d, KMCMC, LMCMC,
                                         # sig2_prior=InverseGamma(11, 5),
                                         sig2_prior=InverseGamma(3, 2),
                                         delta0_prior=deltaz_prior,
                                         delta1_prior=deltaz_prior,
                                         alpha_prior=Gamma(0.1, 10.0),
                                         yQuantiles=[.01, .1, .25],
                                         pBounds=[.4, .8, .05],
                                         similarity_Z=simz)
  s = CytofRepFAM.Model.genInitialState(c, d,
                                        sb_ibp=false,
                                        allow_repeated_Z_columns=allow_repeated_Z_columns)
  t = CytofRepFAM.Model.Tuners(d.y, c.K)
  X = CytofRepFAM.Model.eye(Float64, d.I)
  println("mu*0 truth: $(simdat[:mus][0])")
  println("mu*1 truth: $(simdat[:mus][1])")
  println("sig2 truth: $(simdat[:sig2])")

  _ = CytofRepFAM.Model.compute_marg_loglike(s, c, d, 1.0)
  @time ll = CytofRepFAM.Model.compute_marg_loglike(s, c, d, 1.0)
  println("Loglike with temper = 1.0, marginalized over lam, gam: $(ll)")

  @time ll = CytofRepFAM.Model.compute_marg_loglike(s, c, d, 2.0)
  println("Loglike with temper = 2.0, marginalized over lam, gam: $(ll)")

  return Dict(:d => d, :c => c, :s => s, :t => t, :X => X,
              :simdat => simdat)
end

printstyled("Test fitting repFAM on simulated data ...\n", color=:yellow)
@testset "update W*, r, omega" begin
  config = init_state_const_data(seed=0)
  cfs = CytofRepFAM.Model.ConstantsFS(config[:c])
  dfs = CytofRepFAM.Model.DataFS(config[:d], config[:X])
  sfs = CytofRepFAM.Model.StateFS{Float64}(config[:s], dfs)
  tfs = CytofRepFAM.Model.TunersFS(config[:t], config[:s], config[:X])

  # Do one update for W_star, r, omega.
  println("r init: $(sfs.r)")
  println("W_star init: $(sfs.W_star)")
  println("omega init: $(sfs.omega)")

  CytofRepFAM.Model.update_W_star!(sfs, cfs, dfs, tfs)
  CytofRepFAM.Model.update_r!(sfs, cfs, dfs)
  CytofRepFAM.Model.update_omega!(sfs, cfs, dfs, tfs)

  println("W_star after: $(sfs.W_star)")
  println("W after: $(sfs.theta.W)")
  println("r after: $(sfs.r)")
  println("omega after: $(sfs.omega)")

  CytofRepFAM.Model.printConstants(cfs)

  # Fit model.
  # For algo tests
  # nmcmc = 100
  # nburn = 1200
  # For compile tests
  nmcmc = 3
  nburn = 5
  out = CytofRepFAM.Model.fit_fs!(sfs, cfs, dfs, tuners=tfs,
                             nmcmc=nmcmc, nburn=nburn,
                             printFreq=1, time_updates=true,
                             computeDIC=true, computeLPML=true, Z_thin=1,
                             computedden=true, seed=0)

  outdir = "results/repfam"
  mkpath(outdir)

  BSON.bson("$(outdir)/out_fs.bson", out)
  BSON.bson("$(outdir)/data_fs.bson", Dict(:simdat => config[:simdat]))
end

ext: jl


text:
#RSDeltaSigmaPort: Calculate frequency spectrum information
#-------------------------------------------------------------------------------

#==calcSpecInfo
===============================================================================#
"""
Calculate spectrum information

# Inputs
 - sig: Signal to analyze
 - NTF: Noise transfer function to compute theoretical noise spectrum
 - fband: Specifies (𝑓min, 𝑓max) band limits of the modulator (normalized 𝑓)
 - ftest: test frequency (normalized 𝑓)
 - M: Number of quantizer steps (nlev = M+1)
 - quadrature:
"""
function calcSpecInfo(sig, NTF, fband, ftest::Float64; M::Int=1, quadrature::Bool=false)
	sig = conv2seriesmatrix2D(sig)
	R, N = size(sig)
	if R!=1
		throw("FIXME: replicate hann window rows for R>0")
	end

	iband = round(Int, minimum(fband)*N):round(Int, maximum(fband)*N)
	itest = round(Int, ftest*N)

	Δ = 2
	NBW = 1.5/N
	spec0 = fft(applywnd(sig, ds_hann(N)))/(M*N/4)
	local 𝑓, specdB, SNR, ePSD
	if !quadrature
		𝑓pts = div(N,2)+1 #Only need half the spectrum for display purposes
		𝑓 = range(0, stop=0.5, length=𝑓pts)
		specdB = dbv.(spec0[1:𝑓pts])
#		@show iband .+ 1, itest-minimum(iband)
		SNR = calculateSNR(spec0[iband .+ 1], itest-minimum(iband))

		#Compute Expected power spetral density:
#		Sqq = 4 * evalTF(NTF ,exp.(2π*j*𝑓)).^2 / (3*M^2)
		Snn = abs.(evalTF(NTF, exp.(2π*j*𝑓))).^2 * 2*2/12 * (Δ/M)^2
		ePSD = dbp.(Snn*NBW)
	else
		throw("Not implemented")
	end

	return (freq=𝑓, spec0=spec0, specdB=specdB, ePSD=ePSD,
		SNR=SNR, NBW=NBW, NTF=NTF, M=M, iband=iband, itest=itest
	)
end

function calcSpecInfo(sig, dsm::RealDSM, fband=nothing, ftest=nothing)
	if isnothing(fband); fband = default_fband(dsm); end
	if isnothing(ftest); ftest = default_ftest(dsm); end
	NTF = synthesizeNTF(dsm)
	return calcSpecInfo(sig, NTF, fband, ftest, M=dsm.M, false)
end

#Last line


text:
using LinearAlgebra, ProgressMeter
using TimerOutputs, Sobol

import ParticleInCell: TwoDGrid, landau_sampling!, FDTD, faraday!, ampere_maxwell!
import ParticleInCell: compute_energy, update_fields!
import ParticleInCell.F90: interpolation!, compute_current!, push_x!, push_v!

include("run_function.jl")

t, energy = run( 1 )
@show nstep = 1000
@time t, energy = run( nstep, npm=500 )
show(to)

open("results_f90.dat", "w") do f

    for i in 1:nstep
        println(f, t[i], " ", energy[i])
    end

end

println()


text:
include("../src/ode_tool.jl")
using ODE

tspan=[0,2*pi]
y0=Float64[1 0; 0 0]
function F(t, yy)
    y=yy[1]
    d0=[0 1; -1 0]*y[:,1]
    d1=[1,1]
    {[d0 d1],1}
end
#(T,Y)=ode23_abstract(F, tspan, {y0,0}) ## get samples from the whole ode path
(t1,y1)=ode23_abstract_end(F, tspan, {y0,0}) ## just get the ending value at the terminal time

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
using Test

using CompScienceMeshes
using BEAST
using StaticArrays
using LinearAlgebra

c = 3e8
μ = 4*π*1e-7
ε = 1/(μ*c^2)
f = 5e7
λ = c/f
k = 2*π/λ
ω = k*c
η = sqrt(μ/ε)

a = 1
Γ_orig = CompScienceMeshes.meshcuboid(a,a,a,0.1)
Γ = translate(Γ_orig,SVector(-a/2,-a/2,-a/2))

Φ, Θ = [0.0], range(0,stop=π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]

# This is an electric dipole
# The pre-factor (1/ε) is used to resemble
# (9.18) in Jackson's Classical Electrodynamics
E = (1/ε) * dipolemw3d(location=SVector(0.4,0.2,0),
                       orientation=1e-9.*SVector(0.5,0.5,0),
                       wavenumber=k)

n = BEAST.NormalVector()

𝒆 = (n × E) × n
H = (-1/(im*μ*ω))*curl(E)
𝒉 = (n × H) × n

𝓣 = Maxwell3D.singlelayer(wavenumber=k)

X = raviartthomas(Γ)

T = Matrix(assemble(𝓣,X,X))
e = Vector(assemble(𝒆,X))
j_EFIE = T\e

nf_E_EFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, j_EFIE, X)
nf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, j_EFIE, X) ./ η
ff_E_EFIE = potential(MWFarField3D(wavenumber=k), pts, j_EFIE, X)

accuracy_default_quadrature_nf_e = norm(nf_E_EFIE - E.(pts))/norm(E.(pts))
accuracy_default_quadrature_nf_h = norm(nf_H_EFIE - H.(pts))/norm(H.(pts))
accuracy_default_quadrature_ff_e =
    norm(ff_E_EFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true))

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

    a, b = 0.0, 1.0
    # CommonVertex, CommonEdge, CommonFace rules
    println("Fast quadrule is called")
    tqd = quadpoints(test_local_space, test_charts, (1,2))
    bqd = quadpoints(trial_local_space, trial_charts, (1,2))
    leg = (BEAST._legendre(3,a,b), BEAST._legendre(4,a,b), BEAST._legendre(5,a,b),)

    return (tpoints=tqd, bpoints=bqd, gausslegendre=leg)
end

function fastquaddata(fn::BEAST.Functional, refs, cells)
    println("Fast RHS quadrature")
    return quadpoints(refs, cells, [1])
end

fastT = Matrix(assemble(𝓣,X,X, quaddata=fastquaddata))
faste = Vector(assemble(𝒆,X, quaddata=fastquaddata))
fastj_EFIE = fastT\faste

fastnf_E_EFIE = potential(MWSingleLayerField3D(wavenumber=k), pts, fastj_EFIE, X)
fastnf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(wavenumber=k), pts, fastj_EFIE, X) ./ η
fastff_E_EFIE = potential(MWFarField3D(wavenumber=k), pts, fastj_EFIE, X)

accuracy_fast_quadrature_nf_e = norm(fastnf_E_EFIE - E.(pts))/norm(E.(pts))
accuracy_fast_quadrature_nf_h = norm(fastnf_H_EFIE - H.(pts))/norm(H.(pts))
accuracy_fast_quadrature_ff_e =
    norm(fastff_E_EFIE - E.(pts, isfarfield=true))/norm(E.(pts, isfarfield=true))

@test accuracy_fast_quadrature_nf_e > accuracy_default_quadrature_nf_e
@test accuracy_fast_quadrature_nf_h > accuracy_default_quadrature_nf_h
@test accuracy_fast_quadrature_ff_e > accuracy_default_quadrature_ff_e


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8182327e-1f6a-11eb-30bf-93e0f331ebfe
begin
	using NeXLCore
	using DataFrames, PlutoUI, CSV

	md"""
	## Material Mixture Models

	Compute mass fraction, normalized mass fraction, atomic fraction and other statistics from a mixture of materials including uncertainties.
	"""
end

# ╔═╡ d6c066e0-1f6a-11eb-3cdb-61eb7120583b
md"""
Enter the composition of the material in this text field using the format:

    material1, mass_fraction1[, dmass_fraction1]
    material2, mass_fraction2[, dmass_fraction2]
    ...
    materialn, mass_fractionn[, dmass_fractionn]

where items in [] are optional.
Materials may be entered as chemical formulae like "Al2O3" or as mass fractions like "0.85\*Fe+0.1\*Ni+0.05\*Cr".

$(@bind matstr TextField( (80, 10), default =
"MgO, 0.1933, 0.002\nFeO, 0.0996, 0.002\nSiO2, 0.4535, 0.002\nCaO, 0.1525, 0.002\nAl2O3, 0.0927, 0.002" ))

Material name: $(@bind namstr TextField(default="K412"))
"""

# ╔═╡ 2ed93282-1f6b-11eb-06c2-019bde67147d
begin
	function parse_matstr(nm, str)
		res = Pair{UncertainValues, UncertainValue}[]
		for line in split(str, "\n")
			items = split(line, ",")
			matu, qty, dqty = missing, 0.0, 0.0
			try
				mat = length(items)>=1 ? parse(Material, strip(items[1])) : missing
				matu = mf2comp(mat)
			catch
				matu = missing
			end
			try
				qty = length(items)>=2 ? parse(Float64, strip(items[2])) : 0.0
			catch
				qty = 0.0
			end
			try
			 	dqty = length(items)>=3 ? parse(Float64, strip(items[3])) : 0.0
			catch
				dqty = 0.0
			end
			if (!ismissing(matu)) && ((qty ≠ 0.0) || (dqty ≠ 0.0))
				push!(res, matu=>uv(qty,dqty))
			end
		end
		return length(res) > 0 ? NeXLCore.mixture(nm, res...) : nothing
	end
	nothing
end

# ╔═╡ 3f77a150-1fd2-11eb-03d8-dd33f72d34fb
begin
	pms, mix = nothing, nothing
	try
		global pms = parse_matstr(namstr, matstr)
		global mix=mf2comp(namstr, pms)
	catch e
		print(e)
	end

	if !isnothing(mix)
		md"""
		## Summary
		#### Mass Fractions

		$(asa(DataFrame, extract(mix, MassFractionLabel), false))
		#### Atom Fractions

		$(asa(DataFrame, extract(mix, AtomicFractionLabel), false))
		#### Normalized Mass Fractions

		$(asa(DataFrame, extract(mix, NormMassFractionLabel), false))
		#### Material Statistics

		$(asa(DataFrame, extract(mix, collect(MatStatTypes)), false))
		#### Material Fraction

		$(asa(DataFrame, extract(mix, MaterialFractionLabel), false))

		## Full Covariance
		#### Mass Fractions

		$(asa(DataFrame, extract(mix, MassFractionLabel), true))
		#### Atom Fractions

		$(asa(DataFrame, extract(mix, AtomicFractionLabel), true))
		#### Normalized Mass Fractions

		$(asa(DataFrame, extract(mix, NormMassFractionLabel), true))
		#### Material Statistics

		$(asa(DataFrame, extract(mix, collect(MatStatTypes)), true))
		#### Material Fraction

		$(asa(DataFrame, extract(mix, MaterialFractionLabel), true))

		"""
	else
		md"?"
	end
end

# ╔═╡ 82fd38b0-2039-11eb-25b6-4d2567ec89dd
begin
	md"""
	The full covariance matrix has been written to: $(CSV.write(joinpath(tempdir(),"mixture.csv"),asa(DataFrame,mix)))
	"""
end

# ╔═╡ Cell order:
# ╟─8182327e-1f6a-11eb-30bf-93e0f331ebfe
# ╟─d6c066e0-1f6a-11eb-3cdb-61eb7120583b
# ╟─2ed93282-1f6b-11eb-06c2-019bde67147d
# ╟─3f77a150-1fd2-11eb-03d8-dd33f72d34fb
# ╟─82fd38b0-2039-11eb-25b6-4d2567ec89dd


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QML

using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Circuit
using QuantumCircuits.QCircuits.Circuit: toQiskit
using QuantumCircuits.QCircuits.Qiskit
using QuantumCircuits.QCircuits.Qiskit: qiskit
using QuantumCircuits.Execute

statevectorsim = qiskit.Aer.get_backend("statevector_simulator")

function calcQiskitStateVector(qc)
    qcirc = toQiskit(qc)

    statevector = qiskit.execute(qcirc.qc, statevectorsim).result().get_statevector()
    return statevector.probabilities()
end

################################################################################
#  U3                                                                          #
################################################################################
circ_N = 1
circ_qr = QuantumRegister(circ_N)
circ_cr = ClassicalRegister(circ_N)
circ = QCircuit(circ_qr, circ_cr)
circ.u3(0, 1.3537, 2.5652, 5.5142)

@test sum(abs.(execute(circ) - calcQiskitStateVector(circ))) < 1e-6


text:
using MechanicalSketch
import MechanicalSketch: °, background, sethue, O, PT, finish, empty_figure
import MechanicalSketch:m, color_with_lumin, drawcart, setline, circle, PALETTE
let
BACKCOLOR = color_with_lumin(PALETTE[3], 70)
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_9.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
end
restart()
circle(O, 1.0m, :fill)
sethue(PALETTE[4])
drawcart()
drawcart(p = O + (0m, -3m), α = 30°)
setline(2PT)
drawcart(p = O + (0m, -6m), α = 45°)
sethue(PALETTE[9])
drawcart(p = O + (0m, -9m), α = 60°)
finish()
nothing
end