text:
# Regularized LBM derivation and adaptive time stepping

Hello,

This is my first post on this forum and first of all I want to express my appreciation for this forum, since it has helped me a number of times with several problems related to LBM. For my master’s thesis I am working on the free-surface LBM, which has been popularized by Nils Thurey. Currently I am working on Regularized LBM (RLBM), as explained by Jonas Latt in his thesis [1]. I have two questions about RLBM, one concerns the derivation of RLBM and the other one concerns RLBM and adaptive time stepping.

First the derivation question. For the derivation of RLBM Jonas Latt combines two equations to get to the end result [2]. Since this was not entirely straight forward to me I decided to write out this derivation. This led me to the conclusion that f^{(1)} approx t_i / (2cs^4) * Pi^{neq}_ab and not f^{(1)} = t_i / (2cs^4) * Pi^{neq}_ab as Jonas Latt in his paper states. I have wrote the derivation in Latex and placed it online at http://home.kpn.nl/riji99ou/derivation.pdf (tex-file can be found at http://home.kpn.nl/riji99ou/derivation.tex if you want to add/change anything). Am I making a mistake in my derivation/reasoning?

Secondly, RLBM and adaptive time stepping. According to Jonas Latt’s thesis, we should scale the variables rho, u, and Pi^(1) appropriately when the time step changes. Then we can recompute the distribution functions with the RLBM procedure. As a checking mechanism we can recompute the velocity and density from the rescaled distribution functions, these should be equal to the rescaled velocity and density. However, they’re not in my case, the difference is of order 10^-3. I have posted my code below and perhaps/probably I am doing something wrong. Help would greatly be appreciated.

N.B.: In his thesis Jonas Latt converts u to a dimensionless formulation differently in equation 4.1 and 4.11. I think 4.1 is correct and 4.11 not. Correct me if I am wrong.

[2] J. Latt and B. Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Mathematics and Computers in Simulation, Volume 72, Issues 2-6.

void LBMsim::regularizedChangeTimeStep(float dtNew_p) {
//Determine average density
double rhoTotal = 0.0;
int numFluid = 0;
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID) {
rhoTotal += rho[x][y];
++numFluid;
}
}
}
double rhoAvg = rhoTotal / numFluid;

//Scale global parameters
double st       = dtNew_p / dt_p;
dt_p            = dtNew_p;
gravity_lb     *= (st * st);
double nuNew_lb = nu_p * dtNew_p/(dx_p*dx_p);  //New viscosity
tau             = nuNew_lb / cs2 + 0.5;

//Scale cell parameters
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID || flags[x][y] == IS_INTERFACE) {
double eq[9];
double neq[9];
double piNeq_xx = 0.0, piNeq_yy = 0.0, piNeq_xy = 0.0;
for (int i = 0; i < 9; ++i) {
double fi = f[next][x][y][i];
//Compute equilibrium distribution function
eq[i] = fEq(i, rho[x][y], ux[x][y], uy[x][y]);
//Compute off-equilibrium distribution function
neq[i] = fi - eq[i];
//Compute off-equilibrium momentum flux tensor
piNeq_xx += neq[i] * ex[i] * ex[i];
piNeq_yy += neq[i] * ey[i] * ey[i];
piNeq_xy += neq[i] * ex[i] * ey[i];
}
//Scale off-equilibrium momentum flux tensor
piNeq_xx *= st;
piNeq_yy *= st;
piNeq_xy *= st;
double f1[9];
for (int i = 0; i < 9; ++i) {
//Compute first-order term of distribution function
double Qixx = ex[i] * ex[i] - cs2;
double Qiyy = ey[i] * ey[i] - cs2;
double Qixy = ex[i] * ey[i];
f1[i]       = w[i] / (2.0 * cs2 * cs2) * (Qixx * piNeq_xx + Qiyy * piNeq_yy + 2.0 * Qixy * piNeq_xy);
//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];
}
//Scale velocity and density
double rhoOld = rho[x][y];
ux[x][y]     *= st;
uy[x][y]     *= st;
rho[x][y]     = st * (rhoOld - rhoAvg) + rhoAvg;
}
}
}

//Check
for (int y = 1; y < SIZE - 1; ++y) {
for (int x = 1; x < SIZE - 1; ++x) {
if (flags[x][y] == IS_FLUID) {
double rhoCell = 0.0, uxCell = 0.0, uyCell = 0.0;
for (int i = 0; i < 9; ++i) {
double fi = f[next][x][y][i];
rhoCell += fi; uxCell += ex[i] * fi; uyCell += ey[i] * fi;
}
uxCell /= rhoCell;
uyCell /= rhoCell;
//The following output shows difference of order 10^-3
cout << uxCell - ux[x][y] << " ";
cout << uyCell - uy[x][y] << " ";
cout << rhoCell - rho[x][y] << "\n";
}
}
}
exit(0);
}

Hello,

Thank you for re-reading my thesis carefully. I can confirm the mistakes you are pointing out:

• Equation 4.11 has a typo, and the correct version is in Eq. 4.1
• Whenever there is a relation involving ^(1) and ^(neq) terms, there is no strict equality, but only an approximate equality up to hydrodynamic scales.

As for your code, it seems to me that there is a mistake at the following lines:

//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];

If I understand right, you expect the particle populations at this point to be rescaled, with the new time step. However, although the off-equilibrium population f1 has been rescaled (it is obtained by regularization from the rescaled Pi1), the equilibrium is not rescaled (it depends on velocity and density which have not been rescaled yet). You need to also recompute the equilibrium with the rescaled variables.

I also have doubts about the following line:

rho[x][y]     = st * (rhoOld - rhoAvg) + rhoAvg;

You want to rescale the density like a pressure, right? In this case, I would say that pressure has units of length-square by time-square, and is rescaled as

rho[x][y]     = st * st * (rhoOld - rhoAvg) + rhoAvg;

Hello Bart,

  sorry for bothering you. I have seen you are implementing the free surface LB model proposed by Thurey. Well,

I’m experiencing some trouble with the implementation of this model (the idea would be to implement this model for the palabos library ) and wondering if you could help me by, for example, sharing your code or accepting to take a look at my
test matlab code.

thanks a lot

Andrea

@Andrea
No bother. When I ask others to help me, I should also be willing to help them, so I have cleaned up the code a bit and put it online at http://home.kpn.nl/riji99ou/freeSurfaceLBM.zip In the readme file that is present I have put a small explanation in which files the functionality can be found. The code is written in C++ and is pretty portable, the makefile is written for Linux, and I use OpenScenGraph for the visualization. If you have any questions, feel free to email me (if you find possible errors please email me). If my code won’t help you in locating your error, I’m willing to look at your Matlab code. Also, perhaps we can help each other somewhat.

@Jonas

Actually I still do not quite understand. I read again through your thesis, particularly section 4.2, and wrote out the whole procedure for myself and also came to the conclusion of equation 4.9 that the regularized f^(1) = ti / (2 * cs4) * Qi : Pi^(1). But the part in section 4.2.3 where you say that you can use Pi^(neq) in equation 4.9 as a replacement for Pi^(1) is just something I do not get. It seems like a circle reasoning that will lead to the usual inaccuracies, but the results I obtain in my simulation prove that it’s not. It’s the last part I do not get about regularized LBM.

[quote=jlatt]As for your code, it seems to me that there is a mistake at the following lines:

//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];

[/quote]

Of course, that’s it! Thank you very much! It’s working now. Marvelous!

[quote=jlatt]You want to rescale the density like a pressure, right? In this case, I would say that pressure has units of length-square by time-square, and is rescaled as

rho[x][y]     = st * st * (rhoOld - rhoAvg) + rhoAvg;

[/quote]

Hmmm, you have a point there. Thurey gave the equation in his paper and I didn’t gave it very much thought since it seemed plausible to scale the parts that deviate from the average, so that the average density of the fluid remains equal. Pressure has units of m^-1 * kg * s^-2, so besides scaling with st^2, shouldn’t we also scale with the new mass (which again depends on the density… or am I mixing things up here?). Secondly, scaling it with st^2 gives less accurate results than scaling with st, although doing no scaling at all gives better results than scaling with st, so I am not sure how I should interpret this. I’ll post some results below about how well the mass is conserved after the time step increases with a factor 5, perhaps it will give some insight:

rho[x][y] = st * (rhoOld - rhoAvg) + rhoAvg, gives:

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
TIME STEP FIVE TIMES LARGER NOW
iter=31   dM=-1.44382e-11	  dRho=-1.44382e-11
iter=32   dM=-1.51203e-11	  dRho=-1.68257e-11
iter=33   dM=-1.42109e-11	  dRho=-1.65983e-11

rho[x][y] = st * st * (rhoOld - rhoAvg) + rhoAvg, gives an error st times larger (st=5):

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
TIME STEP FIVE TIMES LARGER NOW
iter=31   dM=-7.38964e-11	  dRho=-7.38964e-11
iter=32   dM=-7.49196e-11	  dRho=-7.61702e-11
iter=33   dM=-7.48059e-11	  dRho=-7.63976e-11

rho[x][y] = rhoOld, i.e. no scaling gives:

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
TIME STEP FIVE TIMES LARGER NOW
iter=31   dM=-1.25056e-12	  dRho=-1.25056e-12
iter=32   dM=-2.84217e-12	  dRho=-4.3201e-12
iter=33   dM=-2.6148e-12	  dRho=-4.09273e-12

Thanks Bart,

         you are too kind...  I will let you have more news soon !

Andrea

Hello,

So, the full argument is as follows. The relation between f_i and Pi used in the regularization formula is exact when you apply it to f_i^(1), but it is approximate when applied to f_i^(neq)=f_i^(1)+f_i^(2)+f_i^(3)+… . Consequently, when you replace the f_i by their regularized counterpart you actually modify their value (and in many cases get a more stable simulation).

I find it intuitive to view this from the perspective of linear algebra in a q-dimensional space (where q is 9 in D2Q9, 19 in D3Q19, etc. I will take q=9 for the sake of simplicity in the following), as it is done in the community of MRT models. f^neq lives in a 6-dimensional sub-space of the original 9-dim space, which is spanned by the vectors Q (3 vectors), h (1 vector) and c (2 vectors), eq. 3.36 in my thesis. During the regularization procedure, the f^neq vector is projected onto the 3-dimensional space of Q, by throwing away the h and c components. I formally identify the result of this with the f^(1) term in the Chapman-Enskog expansion, because the result depends only on a hydrodynamic quantity (the stress tensor), and not on the ghost variables N and J.

You are right, there exist different points of view on this and probably no unique truth. The question is whether the density should be rescaled like, well, a density or like a pressure, which unfortunately is not the same. In grid refinement, this problem is often avoided by using “convective scaling”, i.e. by changing the time step by the same factor as the space step, which results in st=1.

In all other cases, I would be tempted to say that the decision whether to scale like a density or a pressure depends on the physics you want to represent with the numerical model. If the fluid is compressible and you actually want to measure density differences (other than the fluctuation accounting for the dynamic pressure term), then I guess it’s useful to guarantee a continuity of density across a change of scale for space or time. In that case however, and unless I am missing something here, I would say that the density is scale invariant, and remains numerically unchanged when rescaling space or time. Indeed, when I do the series development of the f_i in view of a Chapman-Enskog expansion, I take rho as a grid-independent quantity (“the mass density in a non-dimensional volume”); I am not sure if you can get the Navier-Stokes equations properly if you assume rho to be dependent on dx or dt. This statement seems to be supported by your numerical experiment which yields best results when rho is scale invariant.

Now, if the aim of a simulation is to solve the incompressible Navier-Stokes equations, I have the tendency to view things from a point of view of pure numerical analysis. In a dimensionless version of the Navier-Stokes equations, there is no density: there exists only velocity and pressure. LB is used to solve this equation, and the density is nothing else than a helper variable for solving the pressure equation. It seems to me, then, that it makes sense to scale it like a pressure.

Lattice Boltzmann is a quasi-compressible method for solving an incompressible fluid, and from the “compressible fluid purist” to the “incompressible fluid purist” there exists a whole spectrum of arguments which are valid in their own right, and which probably provide best results in some case or another…

As for the scaling issue with adaptive time stepping, I get it now. And for my problem, I have a compressible fluid, the density therfore remains scale-invariant. Indeed the grid-refinement papers circumvent this problem always by scaling the time step and space step identically. There are by the way only few papers about adaptive time stepping for LBM, as far as I know, only your thesis and Nils Thurey’s paper (and I am pretty sure that Nils Thurey makes a few errors with his equations, but I still have to verify this better and perhaps will post about this issue in the future).

[quote=jlatt]
During the regularization procedure, the f^neq vector is projected onto the 3-dimensional space of Q, by throwing away the h and c components. I formally identify the result of this with the f^(1) term in the Chapman-Enskog expansion, because the result depends only on a hydrodynamic quantity (the stress tensor), and not on the ghost variables N and J.[/quote]

I did get the idea that we want to enforce that the f’s only depend on the hydrodynamic moments (six variables) and not on the total 9 variables of D2Q9. The only thing that I missed is that when you project f^neq to Q that you can identify this with f^(1). Because, as you say it, it only depends on the stress tensor. Perhaps I missed something huge, but I am not aware of the fact that terms that only depend on the stress tensor belong to f^(1).

Secondly, will regularized LBM not get into trouble if we initialize the distribution functions wrongly? That is, the f’s do not comply with f = f^eq + f^(1), and therefore f^neq has errors higher than regularized LBM assumes. Or will this quickly fade away when the simulation runs?

Hello,

If I remember right, this idea has been applied in a paper in which it is called “Mach-number annealing”: A.M.M. Artoli; A.G. Hoekstra and P.M.A. Sloot: Accelerated Lattice BGK method for unsteady simulations through Mach number annealing, International Journal of Modern Physics C, vol. 14, nr 6 pp. 835-847. 2003

The f^(1) obtained by regularization has exact symmetry properties which in the f^neq are only approximate. For example, the relation f_i^(1)=f_opp(i)^(1) (the “bounce-back of the off-equilibrium”) is exact after regularization, which means that the regularized f^neq reproduce these properties of the “theoretical f^(1)” exactly.

To call the after-regularization off-equilibrium populations f^(1) is a formal matter. If they were exactly equal to f^(1) they would exactly match the gradients of the velocity field, which they don’t: this is in any case a numerical approximation. An exact claim however (and that’s what really matters here), is that the variables have been down-projected to a lower-dimensional space, which corresponds to a truncation in the Hermite expansion.

Thank you for all your answers. I think I pretty much get the whole idea of RLBM now and all the terms used in deriving it.

One small last remark about the adaptive time stepping. In the case of simulating a gravity-driven flow, as I do, it might actually not be a bad idea to scale the density with st^2. I am saying this not from an accuracy point of view, but scaling the density with st^2 prevents that the fluid will experience an internal shockwave. That is, gravity causes a pressure difference, which directly affects the density, and if we scale gravity with st^2 we should also scale the density with st^2 to prevent this internal shockwave. This is clearly visible if you visualize the velocity magnitudes of the fluid. Just wanted to add this issue in case somebody might experience some odd behavior with gravity-driven flows and adaptive time stepping.

As for the Mach number annealing paper, I checked it out and it is not really the same as adaptive time stepping. Thanks you the suggestion though.

Now I am working on grid refinement I’ve taken a better look at the momentum-flux tensor (Pi) and how to scale it. In Jonas Latt’s thesis the physical (/dimensionless) representation of the momentum-flux can be obtained from the one in lattice units as follows (Eq. 4.11):
Pi_p = dt_p * Pi_lb

I do not understand this scaling. The SI units of momentum-flux is (http://en.wikipedia.org/wiki/Flux):
N * s * m^-2 * s^-1
With newton defined as kg * m * s^-2 we obtain:
kg * m * s^-2 * s * m^-2 * s^-1
We can simplify this to:
kg * s^-2 * m^-1
Assuming the mass remains constant during scaling, we obtain:
s^-2 * m^-1

With these units we would obtain Pi_p from Pi_lb with the following scaling (which is different from Jonas Latt’s version):
Pi_p = dt_p^2 / dx_p * Pi_lb

I have to make the remark that Jonas Latt in his eqaution deals with the first-order term of Pi (but does that change the units?).

Lastly, to make things more complicated I stumbled upon a paper which uses regularized LBM for grid refinement [3] and it states that:
Pi^(1)_c / (dt_c * tau_c) = Pi^(1)_f / (dt_f * tau_f),
where _c stands for the coarse grid and _f for the fine grid (cells of coarse grid are twice as large as fine grid)

The above equation looks like a third version and I am losing sight here. Can somebody explain how the (first-order term) momentum-flux tensor should be scaled?

[3] Bo Liu and Arzhang Khalili, Lattice Boltzmann model for exterior flows with an annealing preconditioning method, PHYSICAL REVIEW E 79, 066701 (2009)

Anyone?

Besides the question of scaling Pi^(1) I am also now questioning if a possible gravity force should be included in the computation of Pi^(1) if the simulation includes a gravity force. In Guo’s gravity model we, for example, also adjust the velocity by computing as u = 1/rho * sum(ei * fi) + 1/2 * F (instead of the usual u = 1/rho * sum(ei * fi) )

Hello Bart,

you are right, the deviatoric stress tensor Pi^(1) indeed has to be modified if you use Guo’s forcing scheme. In principle, this term is already written down in Guo’s article from 2002 in the second line of eq. (25b) (the term proportional to vF + Fv). This term can be absorbed in the definition of the stress tensor (\sigma) by writing (setting \Delta t = 1): \sigma = - (1 - 1 / (2 \tau)) \Pi^(1) - (1 - 1 / (2 \tau)) (\vec v \vec f + \vec f \vec v) / 2.

Best,
Timm

Hello Timm,

I’m afraid I am not following you. I read Guo’s article more carefully than I first did and equation 25b refers to the resulting analytical momentum equation if Buick and Greated’s model is used (a model that is not entirely correct according to Guo). Secondly, you are mentioning deviatoric stress tensor and I was mentioning the (first-order) momentum flux tensor.
Although I haven’t verified Guo’s reasoning into detail, the following quote makes me believe that the momentum-flux tensor is correctly computed using Guo’s model and his adjusted velocity:
“It is demonstrated that in order
to obtain the correct continuity equation, the ?uid velocity
must be de?ned such that the effect of the external force is
included, and to obtain the correct momentum equation, the
contributions of the force to the momentum ?ux must be
canceled. The method proposed in this paper matches both
conditions, and gives the correct equations of hydrodynam-
ics.”

My question regarding the scaling of the (first-order) momentum-flux tensor is still open. So if anyone has the answer, it would be greatly appreciated. I have done some more thinking about this and Pi^(1) is related to the strain-rate tensor with a constant factor. The strain-rate tensor depends on the gradients of the velocity, i.e. the acceleration, defined in m s^-2. So I would indeed think that Pi^(1) should be scaled as:
Pi^(1)_p = dt_p^2 / dx_p * Pi^(1)_lb

Ah, I made a slight mistake… You are right about the equation.
I try to explain what I mean. The tensor Pi^(1) you are talking about is the second moment of f^(1), as you can see between eqs. (11) and (12) in Guo’s article. In eq. (12), you find a correction term introduced by the force (gravity or any other external force). This equation is the correct result within Guo’s approach.

If you now want to compute the deviatoric stress tensor (gradient of velocity for an incompressible fluid) from it, you have to do it according to $\sigma = - (1 - 1 / (2 \tau)) \Pi^(1) - (1 - 1 / (2 \tau)) (\vec v \vec f + \vec f \vec v) / 2$ (my previous post). This result you obtain when you solve for $\sigma = \rho \nu (\nabla u + (\nabla u)^T)$ and by using eqs. (12), (15), and (16).

Any lattice effects and forces are already included in these equations. By the way: The unit of sigma and Pi^(1) is N/m^2 which is the unit of pressure.

Please apologize if I still have not explained what you are asking for.

Timm

You’re right about eq. 12 for Pi^(1), together with C from eq. 16.

As for the units of Pi^(1), if it’s N * m^-2 and newton defined as kg * m * s^-2 we get:
kg * m * s^-2 * m^-2
Simplifying this to:
kg * s^-2 * m^-1
Assuming mass remains constant, we obtain:
s^-2 * m^-1

So this again yields the same units as I derived above and will also lead to the following scaling to obtain Pi_p from Pi_lb:
Pi_p = dt_p^2 / dx_p * Pi_lb

And as I mentioned in the post above, this is different than what Jonas Latt mentions in his thesis in eq 4.11. Although I think I know something about LBM, I do not consider myself to be more exerienced in this subject than Jonas Latt, which is why I keep doubting about my derivation. And why the work of [3] again has another version of the scaling.

The scaling you are talking about is not correct.
It is the density which has to be constant, not the mass. The reason is that a smaller volume has a smaller mass but the same density. If you follow this approach, you will find that Pi_p = rho * dx_p^2 / dt_p^2 * Pi_lb which becomes, in the diffusive scaling, Pi_p \propto dt_p^-1 Pi_lb.
Now, the confusion should be perfect.

Indeed it’s not getting more insightful, but perhaps after some reasoning we hopefully find the answer.

I think you’re right about that the mass does not remains constant in all cases, although I think this is the case when we are talking about adaptive time stepping, since the volume of the cell doesn’t change then. But you are right that it’s better to talk in the general case. But I’m not following your derivation completely to Pi_p = rho * dx_p^2 / dt_p^2 * Pi_lb and why this leads to Pi_p \propto dt_p^-1 Pi_lb. Could you elaborate this derivation in a few substeps?

Okay, the unit conversion problem… This is how I usually deal with it:
I first write down the unit of the quantity of interest. In this case, it is the pressure or deviatoric stress tensor. The unit is N / m^2. Now, I find the (unique!) combination of the units of density, length, and time for this. So: N / m^2 = kg m / s^2 / m^2 = (kg / m^3) (m^2) (s^-2). I think, this should be very clear.
Now, the next step is to write the pressure in terms of lattice units and physical units, p_P = p_L * C_p. The conversion factor C_p for the pressure then is C_p = C_rho * C_l^2 / C_t^2 where C_rho, C_l, and C_t are the conversion factors for density, length, and time. Not that the conversion factors carry the physical units of the quantity. That’s it. You just have to know the unit of the quantity of interest and the conversion factors for density, length, and time. If the lattice values for density, lattice constant, and time step are all unity (as I always do it), the conversion factor for length is the lattice constant (in meters), for the time is the time step (in seconds), and for the density is the actual fluid density (in kg / m^3).
Is it clear? There is absolutely nothing more about it.

Timm

Hello Timm,

Ah finally I get and I now see that I indeed made a stupid mistake in the derivation in my last post by mixing up Pi_p and Pi_lb. Now how could I do that? Since I just helped someone else with a conversion problem. Anyway, thank you for fixing it and thank you for your derivation, I think I get it now.

So indeed we get Pi_p = Pi_lb * rho * dx_p^2 / dt_p^2

The only issue in my case is, is that I have a gravity force and then rho is not equal to 1 everywhere. As far as I can see it, we then have Pi_p[x][y] = Pi_lb[x][y] * rho[x][y] * dx_p^2 / dt_p^2, with [x][y] denoting the x,y index of a cell (assuming a 2D grid).

I performed a small test with adaptive time stepping of a gravity-driven flow by making the time step five times larger and do the following scaling operations:
rho[x][y] = (newDt/oldDt) * (newDt/oldDt) * (rhoOld - rhoAverage) + rhoAverage; //New rho
//Next compute old Pi^(neq) and scale it appropriately
newPiNeq_xx = oldPiNeq_xx * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt); //New Pi^(neq)
newPiNeq_yy = oldPiNeq_yy * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt);
newPiNeq_xy = oldPiNeq_xy * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt);
//Next compute new distribution functions using the new Pi^(neq)
This yields an average error of 1.1 * 10^(-11)

If we scale Pi^(neq) as follows:
newPiNeq_xx = oldPiNeq_xx * (oldDtoldDt) / (newDtnewDt); //New Pi^(neq)
newPiNeq_yy = oldPiNeq_yy * (oldDtoldDt) / (newDtnewDt);
newPiNeq_xy = oldPiNeq_xy * (oldDtoldDt) / (newDtnewDt);
we then obtain an average error of 4.7 * 10^(-11), which is roughly five times larger. Does this make sense?

As for diffusive scaling, e.g. in grid refinement, we have scale dt_p ~ dx_p^2, and we scale dt_p linearly and dx_p quadratically, then indeed we get Pi_p = Pi_lb * rho * 1/dt_p. Perfect.

Last question for you Timm, the above derivation is different than Jonas Latt’s version (eq 4.11 in his thesis):
Pi_p = dt_p * Pi_lb
and different from the version in [3] that use it for grid-refinement:
Pi^(1)_c / (dt_c * tau_c) = Pi^(1)_f / (dt_f * tau_f),
where _c stands for the coarse grid and _f for the fine grid (cells of coarse grid are twice as large as fine grid)
Was everybody wrong the whole time? Or is there some sense in the above equations?

Okay, in Jonas’ thesis (eq. (4.11)), delta_t is something different, it is not the time step (see eq. (2.75)). This explains the apparently different behavior.

Yes, the density is not exactly constant, but still it is close to being constant. Do not forget: LBM is a slightly compressible fluid solver. There is no one-to-one correspondence here. I usually think in terms of average density when I perform the unit conversion. The fluctuations of the density in the LBM then correspond to the non-constant pressure in the fluid.
Since the pressure scales like the deviatoric stress tensor, the deviations of the local density from the average density should also scale like the deviatoric stress tensor.
In the diffusive scaling, for the pressure, p_P = p_L / dt_P, and the density fluctuations should decrease like dt_P (so that p_P is constant). However, if you keep the spatial resolution fixed, the scaling is p_P = p_L / dt_P^2, and the fluctuations should vanish quadratically. However, the spatial error is not reduced at the same time, and this will not work infinitely. You also have to change the relaxation parameter which will lead to additional problems.

The errors for the stress tensor, are they absolute or relative? Please define them correctly.

What do you mean by reference [3]?

Timm


text:
# Math Help - Mod 15

1. ## Mod 15

Q: If 5x=5y (mod 15), then could someone please explain to me why it's NOT true to say that x=y (mod 15) but it IS true to say x=y (mod 3)? Just a little confused as to why I can't divide both sides by 5 in the first instance to give x=y (mod 15) and yet I can for x=y (mod 3).

2. Would a counterexample be enough? Try x=4, y=7.

3. Hello,

$5x=5y \bmod 15 \implies 5x-5y=15k \quad k \in \mathbb{Z}$

Simplifying : $x-y=3k$, that is to say $x=y \bmod 3$

This is partly because if you assume that x-y=15k' (which would yield x=y mod 15), it means that k was a multiple of 5, and nothing lets you say that.

format: plain


text:
Confession of a Haskell Hacker (r6.ca) 163 points by malloc47 on July 8, 2012 | hide | past | web | favorite | 124 comments

 Your description on deleting files not being protected by the type system is the result of side effects at runtime.The point of the OP was discussing the safety that can be achieved from pure code.
 I feel like I addressed this fully in my comment. I agree with the OP - you can do great things with the Haskell type system, if you're writing pure code and if you're writing sufficiently polymorphic code. I think this is awesome!However, aside from a few libraries and a smattering of very heavily theoretical work, no code satisfies both of those predicates. If you want to write a 'real world' application in Haskell, you're going to have to get your hands dirty with IO and with concrete types, which means that the ability of the type system to prevent you from shooting yourself in the foot is severely curtailed.
 I find this happens a lot with Haskell. I can write hundreds of lines of code and have them work perfectly as soon as I get a error-free and warning-free compilation.The type system is definitely a big part of this, but almost as important are algebraic data types and compiler checks for exhaustivity when scrutinising values of such types.
 This just screams "bad development practice", IMO.Maybe I've only worked with amazing super hero developers (not likely), but the type of bugs avoided by Haskell's type system just aren't a big problem in any of the projects I've ever worked on.The difficult bugs are almost always mistakes in the requirements and design errors, which become even more difficult to fix when you have a bunch of code.On the other hand, I guess I could see type errors being a problem if you often write hundreds of lines of code before trying to compile and run it...
 Bad requirements and design errors can't be fixed by a perfect programming language.Those are certainly the biggest problems in building software, but not the types of problems this particular article is talking about overcoming.(I think we're both in agreement on this)
 This is supported by the parametricity theorem. That is a big word, but it boils down to this: let's say you wrote an identity function of type "a -> a", and it passed the type checker. Then it is correct: you simply can't do much with a value of type "a", because you don't know anything about it.If an identity function is too simple, consider "compose :: (a -> b) -> (c -> a) -> (c -> b)". I think it can be proven that if you write an implementation that passes the type checker for this signature, the implementation is necessarily correct.
 That's not quite true, it must either be correct or bottom.I can write: id :: a -> a id = id
 This is why Haskell is evil and you should use ML.
 This is why Turing-complete languages are evil and you should use Agda.
 Yes in ML you'd have to eta-expand first: val id : 'a -> 'a fun id x = id x
 You can still write x=_|_ in ML, you just have to be subtler about it.
  id x = id x
 For certain definitions of correct... id :: a -> a id a = case (a,unsafePerformIO (print "rm -rf *")) of (x,()) -> a
 Note the "Safety Inferred" seal of approval on the upper right hand corner of the hackage package http://hackage.haskell.org/packages/archive/lens-family-core... This means it survives a recursive check of dependencies that precludes that kind of nonsense.
 I hadn't thought of that. Types serve as a good safety net, but they're not foolproof. Expecting something like that to be correct because it passed through the type checker is like expecting to be served at a restaurant dressed in a gorilla suit because you got through the door.
 And of course it is fun to think about total functions that are "effectively" bottom... id :: a -> a id a = case (a,ack 4 3) of (x,0) -> x (x,_) -> x where ack 0 n = n+1 ack m 0 = ack (m-1) 1 ack m n = ack (m-1) (ack m (n-1))
 Hi. I make a living writing code. Also, I'm bad at math other than bit arithmetic and calculating tips. (Edit: I guess I sort of understand sets and basic set operations)What in the hell are you two talking about and why shouldn't I take this as a sign that I should avoid Haskell at all costs?
 You can prove properties in Haskell using the type system, (provided that you can show that the expression you are proving properties about actually terminates).I was merely pointed out that last part in the most succinct way I knew how.You should definitely not avoid Haskell at all costs. I would characteristic programming as the process of constructing systems and managing complexity. If you have a greater exposure to languages and paradigms then you have more tools at you disposal for doing your job. A programmer may be able to get by without knowing anything about Haskell or functional programming in general, but that's only because they don't realize when they are kludging something that is so simple if only it was thought about differently.
 I agree that laziness turned out to be a bad default for the reasons you mentioned (other Haskellers may disagree with me on that point), but I don't think that means that it was a mistake to pursue it. It forced the designers to build in purity which turned out to be immensely valuable. Having laziness as a default is also not the kind of decision you can easily revert.You're right to say that the field of static analysis has been very active. What you're forgetting is that this has been in no small part due to type theory.I think the reason you don't find people talking about actual code, as you put it, is because actual code and so called mathematical gibberish are one and the same thing.Lisp is based on the idea that code are data are equivalent, Haskell is based on the idea that code and mathematical gibberish are equivalent.I didn't intend to imply that you didn't know the benefits of functional programming. I could have worded that better. Do you know the benefits of typed functional programming?Having recently delved into lisp, I do miss the static analysis that goes on before the program starts. I realized that I used types to plan out the program before I started coding and I miss having the type-checker picking holes in my plans. I'm determined to stick at it though so that I can see what I'm missing in terms of meta-programming and other dynamic features.
 >Do you know the benefits of typed functional programming?Yes, that's how I know what ML and the tragedies of tagged pointers are.SBCL supports type annotations, and Racket has a Typed Racket dialect that does type checking insofar as you've provided the types.Part of the problem that Lispers and Schemers have is that they can add the static analysis to their programs anytime they want through the macro system (hence typed racket). They don't have to hack up the compiler or anything.The other problem is cultural. Lispers resent implementation-level obligation to a set of rules that were made ignorant of any given problem they may be solving. What if you need to resolve something at runtime? What if you need to resolve something at compile-time? Lisp lets you choose where each should happen.I feel similarly about purity. I think it makes for a nice default, but I think the degree to which immutability is forced on you in Clojure and Haskell are pretty damned awkward for expressing a variety of problems and not just ones related to IO. It seems especially silly since a lot of code in Common Lisp will by default be pure because it'll just return results based on inputs absent mutation.Regarding your delvings into Lisp, since you're not unfamiliar with the benefits of functional programming I'd like to recommend a book that will help you break ground with a dynamically typed, multiple-dispatch Lisp with an object system. Get "Art of the Metaobject Protocol". Hopelessly specific to Lisp but glorious nonetheless.The book doesn't get as deep as I'd like, but it's considered one of the best books on the subject nevertheless.
 One of the things I don't like about Haskell is that it has a really complicated type system because it isn't (and couldn't realistically become) dependently typed. This is one of the things that drew me to lisp, you can choose which type system you use.I must disagree with you on purity though. Purity allows you to apply equational reasoning, it allows compilers to apply optimizations that wouldn't otherwise be possible. It makes multi-threading, if not possible then substantially easier. You do get some of these advantages from having purity by default, but if it's not enforced, you can't get most of them. For me, I see this as an acceptable compromise (actually I don't see it as a compromise at all, if all I got from purity was equational reasoning I'd still use it), but as you point out, the culture amongst lispers is different.Thanks for the book recommendation, I must admit that I already have a number of books to get through before I'll be able to look into it (including LiSP and "Compiling with continuations"), but I surely will do.
 My dream programming language is actually a Lisp with static-by-default typing that was powerful and useful. If you end up implementing something, PLEASE post it on HN or let me know, I'd love to check it out. :)This is how I know about Typed Racket despite not being that keen on Scheme dialects (other than call/cc), I was very excited to see how compile-time typing could work in a Lisp/Scheme.
 > If you don't know the c2 page for this, you should feel deep shame.Consider me shamed, mind sharing a link? I know the quote from Alan Perlis you're talking about, but the Google search results for "perlis languages" is dominated by Fogus's blog post with the same title, and your comment is the first result for "perlis languages c2", followed by a few other wiki pages.
 Types are propositions, and programs are proofs. If your program satisfies the type, then it is a proof of the proposition.
 That is true, but the proposition is one that is encoded in the type. There are many derivations of certain proofs. Consider: map :: (a -> b) -> [a] -> [b]  And [a] is shorthand for a recursive type: \forall a . \mu l . (a * l) + ().And the translation of types to propositions means that the sum type -- either I have () or I have a list, translates to a logical connective of "or". The product type -- I have an a and the rest of the list -- translates to "and" .The proposition is therefore that, for all a, for all b, if a implies b, then true or (a and list of a) implies true or (b and list of b).There are several proofs of this proposition, and some of them are useless.The standard one: map f [] = [] map f (x:xs) = f x : map f xs  Is what we want.There's another one, though: map f [] = [] map f (x:xs) = [f x]  And many others that aren't so useful as the standard one.Anyway, this has been a long and roundabout way of saying that, while a good type system give you a lot of nice static properties about a program, you simply cannot ignore dynamic semantics and rely on the types. Well-typed programs are correct in that they will not reach a state where evaluation is undefined / get stuck ("I'm supposed to do what? These are Strings, not Ints!"), but they do not imply correctness of the program, even though they do provide a proof of a proposition!
 In Haskell all types are inhabited, so every proposition has a proof. This makes Haskell a pretty useless logic.
 Some of them are only inhabited by bottom, no? I cannot write down a function f :: a -> b besides: f :: a -> b f = undefined
 There are a lot of other Haskell features you can use. Term-level recursion works: f :: a -> b f = f  Type-level recursion works, even without explicit term-level recursion: data T a b = C (T a b -> (a -> b)) q :: T a b -> (T a b -> (a -> b)) q (C f) = f w :: T a b -> (a -> b) w x = (q x) x f :: a -> b f = w (C w)  This is basically encoding Russell's paradox at the type level. You can write out f explicitly, just so that it doesn't look like you might somehow still be applying w to itself recursively: g :: a -> b g = (\x -> ((\(C f) -> f) x) x) (C (\x -> ((\(C f) -> f) x) x))  There are even ways of doing this using highly impredicative definitions involving GADTs and type families, without involving any explicit recursion at all:http://okmij.org/ftp/Haskell/impredicativity-bites.htmlThe GADT example no longer seems to work, so changes to GADT type checking might have eliminated it, but the type family one still does.
 It doesn't affect your objection to MaleKitten but not all types are occupied in Haskell. You are thinking of kind * By contrast, {-#LANGUAGE DataKinds#-} a :: 'Nothing a = a  and {-#LANGUAGE MagicHash#-} import GHC.Prim a :: Int# a = a  don't typecheck, though I cant say I understand the ins and outs of the latter.
 Find the type of:fix x = x x
 From the article: "The types are so polymorphic that I conjecture that there is only one way to write functions matching the required types such that all parameters are used non-trivially and recursion is not used."He is not arguing "it works if it compiles" in general. He is arguing "it works if it compiles", because types are so polymorphic. Polymorphic part is actually important. If you don't get that part, you are missing the point.
 Ah, a Haskell circle-jerk waiting to happen and a confession that's really a way to act superior about Haskell's type system. Listen, it's a great language and I've spent a few years deep in it, but I find this post arrogant and not even accurate.The libraries he is talking about are for data types, something that the type system (almost by definition) is perfectly adequate for testing. So yes, he's tested it by compiling it. But for any real world Haskell program the type system alone is not enough.
 Wow, you are quick to accuse the community of being a "circle-jerk" even before most people commented. Looking through the other comments on this thread, most are either against the practice or qualify it. And the qualification is completely reasonable--given sufficiently general types, the space of possible implementations becomes very small which means you are increasingly likely to be correct just if your code typechecks.The article never claims you can always use the type system in place of tests. It doesn't even claim that it is a reasonable practice. The only thing it says is that the author did rely on the type system and that he thinks it is reasonable given sufficiently polymorphic types.In short, you're knocking down a straw man and using it to accuse the entire community of being an arroganet circle-jerk.
 I think your reply is inappropriately negative ("circle-jerk", "real world"), and that the negatively of your post detracts from what value might be derived from the article.The author of the post provided an amusing example of a case where he posits the type system is sufficient and valuable. He didn't claim that the type system would be sufficient for a complete real-world application.
 Well, I guess I'm just old enough and ugly enough to call things like I see them. Haskell is a neat language, but it is frustrating to listen to these sorts of brag-posts about the language because they're filled with dangerous hubris. The type system is amazing for certain kinds of problems. Hey, if you're writing a compiler then Haskell is probably the best language out there. Or in memory data structures. Hm. Or parsers. . . . It's a great language for CS research.The term "real world" is not meant to slur, it's meant to be accurate. I've shipped code written in Haskell and I deeply understand the huge pain in the arse it is to use in production for applications that are not toys.
 Sorry to burst your bubble, but I wonder if your experience is clouded by the (past) inadequacy of your skill set and learning level. We use Haskell in production across a number of domains and it has produced solid results. So do a number of other very serious companies out there (as an example, most banks have varying degrees of Haskell in production for critical systems).Also, how are parsers or in-memory data structures NOT real-world tasks? You never had to parse some bit-stream or stage a sophisticated-enough computation pipeline that required lots of different, intermediate data types?I'm not sure why you're seemingly dying to make negative blanket statements based on how "it was a pain in the arse" for YOU.
 "Dangerous hubris" is an extremely condescending phrase. Do you think the author does not realize that the type system is only sufficient in certain, very specific cases? That's exactly what he mentions at the very end--the reason he trusted the type system so much is only because he believes there is only one non-trivial, non-recursive implementation with the requisite types. He qualifies exactly what he means and, moreover, does not even defend his behavior. And that is "dangerous hubris"?Also, there is certainly plenty of real-world code I would not like to write with Haskell--drivers, for example. However, it seems great for web-apps, which are very "real world" and, naturally, very important on HN.In my personal experience, it's also surprisingly good for quick command-line scripts in place of Bash or Perl. Writing small scripts for automation and command-line clients for my existing tools was as easy or easier in Haskell than in any other language I've used.
 I don't think you realize just how condescending and full of arrogance your posts are.You are implying that computer science is both useless and only consists of implementing toys. What makes your problems "real world" and the problems that haskell is good for not "real world"? What makes something "not a toy"? Parsers and compilers are toys?If you want to start a discussion about things that make haskell difficult to work with on a daily basis, please do so. Otherwise, kindly stop throwing around insults.
 Can you share some examples where you've found it to be a huge pain in the arse to use in production?
 If you're using the term “real world”, it is extremely unlikely---like sha1 hash collision unlikely---that you are being accurate. In all cases where I can remember its being used, it's being used as a way to avoid a concrete and accurate description of fact, and to add a touch of disdain for something else.
 >I've shipped code written in Haskell and I deeply understand the huge pain in the arse it is to use in production for applications that are not toys.And yet plenty of people (myself included) are using haskell for real world, production, non-toy applications all the time. So either all of us people using it are mentally ill and just imagining up non-existent production haskell code, or your experience reflects your ability, not the language's capabilities.
 > not even accurate.What is not accurate in that post?
 I explained that in the second paragraph.
 The author didn't claim that the system is alone sufficient "for any real world Haskell program", so your refutation appears rather spurious.
 I explained that he _did_ test the program because he was writing a type. The type system by definition checks types.
 I see what you mean, but it's rather bizzare to classify type checking as a form of testing, even if type system is very strong. Types and tests are different ways of specifying behaviour: tests are existential guarantees (output of the code sometimes has some properties), types are universal guarantees (output of the code always has some properties).In an ideal world, tests would be subsumed by types, but in reality tests can cover things that types do not.
 > In an ideal world, tests would be subsumed by typesIt would be literally impossible for this to happen in a Turing-complete language, though, because a hypothetical 'perfect' type system would be able to solve the halting problem.The Haskell 'Bottom' type alludes to this in its documentation, actually: http://www.haskell.org/haskellwiki/Bottom(Essentially, Bottom is a member of every type, or else the compiler would need to distinguish between a halting program and a non-halting program before accepting the input program).
 Some people argue that in an ideal world we wouldn't be using Turing-complete languages. Nothing like a total functional programming language to make life better :).
 Some people argue that totality isn't enough, you really want a feasible (polytime) language.http://www.di.unito.it/˜gaboardi/papers/BaillotGaboardiMogbi...;-)
 I thought about total functional programming, like Agda. But even if you have Turing completeness, you can add time constraints to tests and everything will be decidable.
 I think the interesting point here is that Haskell provides such a high static assurance out of the box that what you write is correct that this can happen. You'd never hear of a ruby programmer releasing a gem without trying it.
 With the serious caveat that the program you are writing must only be vulnerable to the kind of errors which Hindley-Milner can show aren't present. ie: no IO, no non-deterministic concurrency, no non-total functions, etc, etc.
 I wish this caveat were attached every time someone mentioned the "once it compiled it worked first time" bit. Not because there's no truth in the works-if-compiles idea, but because getting smug about it is sure to bite you in the ass sooner rather than later.Now that I've said that, I'll state what should be obvious (and probably is to many/most): failure to compile due to type errors is a big clue I am thinking about something incorrectly or incompletely. It's not like I just twiddle stuff until it compiles. I look at it and say, "Duh! That was dumb of me!" It's like using Unit Analysis working out a long physics problem... errors in your thinking about it are very likely to show up in unit analysis. The chance of you coming up with a bogus equation that passes unit analysis is pretty small. When it doesn't pass, you have a good look and find where you went wrong rather than merely tweak.
 That little connection sent me on a Google quest.Interesting paper along those lines:"Types for Units of Measure: Theory and Practice": http://research.microsoft.com/en-us/um/people/akenn/units/CE...
 Oh! I am so glad you looked that up. Thank you!
 You forgot the bug: "doesn't do what you want, but does it perfectly".
 In addition to that, it's often difficult to predict the runtime properties of a Haskell program. Due to it's laziness, one can often run into space leaks where thunks are not evaluated (yet):
 That's still a significant portion of the code-base though. In almost every other language you'd have to write tests to accomplish the same thing. It's certainly possible to compose functions that return something of the correct type but not the correct value (in fact, that's trivial to do). Those should be easier to spot/debug/test than the errors the type system checks for you.
 it's a matter of how much you put into the type system: add :: Integer -> Integer -> Integer add x y = x * y  in any moderately complicated code, you need to test. types only go so far.
 To supplement "it's a matter of how much you put into the type system":This error can be caught if you put units of measure into type system. F# can do this, and addition function will have type int -> int -> int while multiplication int -> int -> int; for example, int -> int -> int might be a function taking velocity and time, and returing distance. I think Haskell has libraries for this.Similarily, add x y = x - y would be detected if (-) required a group, while (+) worked for monoids only.However, it would be rather hard to create a type system detecting a typo in:quadruple x = 3 * x
 > it's a matter of how much you put into the type system:It's also a matter of how much polymorphism you have. There are lots of functions Integer -> Integer; there is only one total function a -> a. (This is the parametricity theorem mentioned by sanxiyn in http://news.ycombinator.com/item?id=4215055 .)
 There are lots of different functions a -> a. One could take a few cycles to complete, whereas others could take a day or week to complete. One could use a few bytes of space, whereas others could use megabytes or gigabytes.
 I think that this is a different meaning of 'different' from the one usually recognised in computer science; it's my understanding that what usually matters is extensional equality (are the same outputs (eventually) produced for the same inputs?) rather than intensional equality (are those outputs produced in the same way?). Otherwise, all compiler transformations would be illegal.(You could argue that a (correct) transformation that sends space usage way up may actually change even extension equality, in the sense that it could result in the transformed algorithm failing on a given piece of hardware when the original algorithm would have succeeded. However, my understanding is that the state of 'realistic' computer science, where hardware limitations are part of the reasoning, is in a much more primitive state than 'idealised' computer science.)
 Opaque types like "Integer" don't give a lot of assurance, so either you get a stronger type system (Agda-style) or you have to test.But complicated code may work with non-opaque types, and so it could still be reliable without tests.
 I have some open source projects that I've released on multiple occasions without having tested them. I think people make the assumption that things that have adequate documentation and comments/support tickets, etc. that indicate use, that they are tested before release. That is simply untrue. You get what you get. That is true in the world of both paid and unpaid open-source and closed-source software and well as life in-general.Something that helps in this regard is travis-ci. It's a free CI server and if you have at least some level of testing, you have some level of confidence.
 Why did you write it if you have never run it? Was it just for people who read your article and might want to run the code?
 Probably as an intellectual exercise.
 But what point is there to an intellectual exercise if you don't know whether you've done it right in the end?
 "Beware of bugs in the above code; I have only proved it correct, not tried it." -- Knuth
 That's the point of the article. The claim is that he knows he has done it right because the declarative types are simple enough to understand as a whole, and they enforce and ensure correctness, verified through compilation.
 Until the code has actually been tested, he doesn't know that it works, he just thinks that it does. The above point is that he has not completed the exercise, by verifying what he thinks.
 > Until the code has actually been tested, he doesn't know that it works, he just thinks that it does.After the code has been tested, he still doesn't know that it works, only that it sometimes produces the expected output. (As a mathematician, it's a lot easier for me to trust a proof than an empirical verification; but I agree that what the HM type system proves is not always what one wants it to prove.)
 I agree. I claim that he knows that it works under the circumstance that he has tested, and can be more confident that it will work in circumstances similar to what he has tested.
 True. But with very general code, like this library, it would be difficult for unit tests to cover more than a small fraction of the domain.
 The code was tested by the compiler, reducing the failure space to whether or not his declarative types are an adequate Proof.It's an amusing anecdote. Even when using Haskell, there are some obvious constraints to relying on the compiler as the sole validation of correctness.
 That's not good enough; you're still testing a proxy behavior, not the behavior under question. The hypothesis is, "Haskell's type system is robust enough that one can code an application through type checking alone." Doing the hypothesis is not testing the hypothesis; assuming the hypothesis is correct to demonstrate that the hypothesis is correct is begging the question.In particular, the kinds of errors that still may exist are ones where the implementer has a fundamental misunderstanding of what needs to happen. The code may correctly implement what he wanted it to do, but what he wanted it to do is wrong.If you still object, consider: I'm describing the scientific method.
 Consider that his Haskell program (ostensibly) relies on mathematical proofs, a chain of them. The compiler consumes and validates this chain of proofs and checks whether or not the program is "true." It's hard to imagine why you'd use the scientific method to verify that 1 + 1 + 1 == 3; basic math and logic--- something computers are good at--- are sufficient.As he and others have pointed out, the author was pretty up-front about the scale (relatively small; looks like a few hundred LOC) and nature (implementing an extant, well-defined API) of the problem. It's not an application as you suggest, and the author doesn't suggest it's a good idea for applications, either.The conclusion isn't "use Haskell & never test your code!". After all, he explicitly offers no conclusion. :) It's just an example of how, sometimes, you can be reasonably certain that your Haskell program is correct when it compiles.Myself, I'd still run it, or at least QuickCheck it. I don't trust myself not to write logically sound, completely incorrect code.
 As an exercise, I think it's great. But my point is that he has not completed the exercise.
 Once you've figured out on paper that 1 + 1 = 2, you don't necessarily have to put one apple next to one apple and count two apples to “complete the exercise.”
 No, you don't. But software is not that simple.
 I object because you're attempting to extrapolate an axiom from an anecdote merely intended to demonstrates the value and power of the type system.
 My point is we have not yet verified if the assumed conclusion of the anecdote is correct.
 I actually encounter this a lot with Haskell. Once you taste the power, it isn't good enough. You want more. Dependent types, etc... But the truth is, it is still a large step up from where we were before in my opinion.
 The compiler can only test what he writes - not what he intends to accomplish.
 Ah yes, the old question of "what is knowledge?".
 Well, unless he has a proof that the code works. I think the idea is that, in this case, the types are so polymorphic that they do constitute a proof.
 The insight behind Knuth's glib "I have only proved the program correct" quote is that proofs can still contain our untested assumptions. It's possible for someone to carry over an untested and wrong assumption from the implementation to the proof without ever finding out it is incorrect.
 Of course, the same is true for tests. If you prove the wrong thing, your code obviously won't be correct. If you test the wrong thing, your code obviously won't be correct.However, if you prove the right thing, your code is correct. If you test the right thing, you code is correct in the cases you tested.And yet, for whatever reason, only having tests is fine but only having proof isn't.
 There's no way to prove that you proved the right thing, unfortunately. And, of course, your proof may not be correct. (Proof checkers can help, and I would have confidence in the correctness of a type-check as proof.)Only having tests is "fine" for most circumstances because, unfortunately, proofs are intractable for most systems. During a software engineering talk I recently attended, an interview candidate stated that the largest proved-correct code base is on the order of 10,000 lines of code. Most people don't have the resources to do even that.But, philosophically, I always want tests for the same reason that I run experiments: I want empirical evidence, not just reasoning. Empirical evidence increases my confidence that we have accounted for what we think we have.
 That's not very surprising; Haskell's type system is really that ridiculously powerful. But I suspect many hackers working with other languages would be able to state the equivalent: "I released a library which I only unit-tested, without writing any helper project that actually uses it".
 except that they had to write unit tests - whereas the compiler/type system automatically checked the type safety of the code this guy wrote.(Not that writing unit tests for haskell code is a bad idea - QuickCheck is WAY more powerful than junits, for example)
 QuickCheck and JUnit have different purposes. QuickCheck is not a unit testing framework.
 > QuickCheck and JUnit have different purposes. QuickCheck is not a unit testing framework.Indeed. What I meant was that basically: P(HaskellCodeCorrect|Compile+QC+UnitTest) =~= P(OtherCodeCorrect|Compile+UnitTest)  while the distribution of effort between compilation (if any) and unit tests in other languages are quite different than in Haskell, i.e. skewed towards tests in the former and getting code to compile in the latter case.
 And that's something he's proud of?If I never ran my code, it wouldn't have any bugs either, no matter what language I was using.I guess it does help explain all the half-assed packages on Hackage.
 Hindley-Milner is powerful, but there are usually multiple ways to construct a value of the correct type. Only one of those ways is the correct way, and I'd want to test my functional code to make sure that this was the way that I chose.
 As sanxiyn points out in http://news.ycombinator.com/item?id=4215055, although there are certainly many ways to produce a value of any concrete type, sufficient polymorphism will often guarantee (via the parametricity theorem) that there is only one total function of a specified type.
 nothacker, you are "dead"/"ghosted". I am reposting his comment here:nothacker 42 minutes ago | link [dead]I have some open source projects that I've released on multiple occasions without having tested them. I think people make the assumption that things that have adequate documentation and comments/support tickets, etc. that indicate use, that they are tested before release. That is simply untrue. You get what you get. That is true in the world of both paid and unpaid open-source and closed-source software and well as life in-general.Something that helps in this regard is travis-ci. It's a free CI server and if you have at least some level of testing, you have some level of confidence.
 Minor aside, is "I conjecture" valid English?I'm not criticising, just curious and in the pub.
 'Conjecture' is a verb as well as a noun. But I'll admit it sounds weird — by French inflection one would expect "I conject", which seems to be obsolete...
 Yep. Dictionary.com calls it obsolete. Webster wants me to pay, and OED can't find it.conject late 14c., obs. verb replaced by conjecture (v.). Also in form congette.
 my concise oxford dictionary (dead tree, c 1984) lists conjecture as a verb (not marked as obsolete).
 It's the word "conject" that's obsolete, not "conjecture".
 oh, sorry, mis-followed the conversation.
 I believe so.
 Very common in mathematical circles. Actually, the entire first page of Google's results for "we conjecture that" appears to be math.
 As someone who is interesting in writing more Haskell, I've heard all the "if I get it to compile it works perfectly the first time", and I wonder if people are missing out on what proper TDD can bring them.At first TDD becomes a way to assert your code does what you expect. After a few years you begin to use it as a design technique, and the correctness argument becomes less and less of a reason for using it. It's more of an (albeit nice) side effect at that point. I can write code that I test after the fact and still get the same correctness benefits, but the feedback I get from testing my design is gone.Maybe I'm missing something, but I don't know if just having an excellent type checker is enough to provide the same quality feedback loop as well executed TDD does.
 You can use the types the same way--define your types and write type signatures for your functions before implementing them.I've found this particularly useful when writing really confusing code (like a Prolog interpreter, for example). Coming up with the types on a function clarifies what it has to do and shrinks the space of possible solutions significantly.Types are a symbolic way to reason about code. Using types, you can get additional insight into how your function can be written just by following some simple rules. A good example of this is realizing that a function you're working on is actually a specific version of a more general function--you can see this if the type signatures look similar. This lets you reuse very high-level, generic code easily.Going back to the Prolog example, I was having some trouble figuring out how the resolution algorithm should work. Then I realized that the sub-part I was having problems with was just special type of fold. This helped me get a very concise version of the function by reusing an existing library function.So types can provide exactly the same sort of benefits as TDD.
 This is why Hackage is full of junk and Haskell is unapproachable for real-world projects. Everyone's throwaway toy file is published as a package, and the useful well-documented packages get lost in the mix.
 Ignoring the baseless FUD in the first part of your post, your last point is actually valid: it is hard to find good, maintained packages on Hackage. Right now the standard solution seems to be just asking somebody with more experience, but this is clearly not scalable.Happily, a new version of Hackage creatively called Hackage2 is being worked on which should ameliorate this problem (along with some other improvements). Improvement is just around the corner.
 Hackage2 has been around the corner for years. The "public testing instance" doesn't exist. http://hackage.haskell.org/trac/hackage/wiki/HackageDB/2.0If someone's actually working on it again, that would be delightful.
 As somebody who knows Haskell only superficially, this line of code and the author's statement makes me wonder if people can read their own or somebody else's code after some time.
 Yes, you can. I'm not a Haskell guru, and I haven't played with it in some time. I can read that line of code just fine. I don't actually know what it's doing, because I don't know the functions involved, etc., but the form of the line is easy enough to grok. Like any other language, it looks mysterious if you haven't read and written enough code in it, or in a language with a similar style.
 "I don't actually know what it's doing"Sorry, but that means you couldn't read it.
 In a very real sense, that's true. But...I don't know what it's doing but I know how it's doing it, so to speak. With languages in which I'm quite fluent (Haskell not so much), I can sometimes spot problems in the code of someone coming to me for help without knowing anything about the libraries, et al, because even though I "couldn't read it" in your sense I could read it in the other sense.
 No it doesn't: foo(bar(x))  Can you read that code? Yes? But you don't know what it is doing, because you don't know what the functions foo and bar do. Same deal with the haskell code in question. Just because you don't know what Compose and getCompose do, doesn't mean you can't read the code.
 I don't feel it's so clear cut. Haskell fuzzies the difference between syntactic and semantic elements. What in your piece of code is the parentheses might very well be defined as functions in Haskell. For instance, the ($) and (.) functions, but also I would call StateT and perhaps Compose and WrapMonad (while not knowing what they do) such things too.So, I still feel that understanding a piece of Haskell means that you will need to understand the individual functions at least superficially.  You feel wrong. There is no such fuzzy. StateT, Compose and WrapMonad are functions. They are just like foo and bar in my example. Just as in a C-style language you need to know the parens mean function application, in haskell you need to know spaces and ($) are function application, and (.) is function composition. The code is perfectly readable even if you don't know what the functions being used do.
 Makes perfect sense to me.

Applications are open for YC Summer 2019

Search:


text:
# Fortran/language extensions

Like several other languages, Fortran 90 and newer supports the ability to select the appropriate routine from a list of routines based on the arguments passed. This selection is done at compile time and is thus unencumbered by run-time performance penalties. This feature is accessed by use of modules and the interface block.

In the following example, a module is specified with contains a single logical function which can handle arguments of various types.

module extension_mod
implicit none
private

interface f
module procedure f_i
module procedure f_r
module procedure f_z
end interface

public::f

contains

function f_i(x) result(y)
integer,intent(in)::x
integer::y

y = x**2-1
end function f_i

function f_r(x) result(y)
real,intent(in)::x
real::y

y = x**2-1.0
end function f_r

function f_z(x) result(y)
complex,intent(in)::x
complex::y

y = x**2-1.0
end function f_z

end module extension_mod

A program which uses this module now has access to a single logical function f which accepts arguments which are of the integer, real, or complex data types. The return type of the function is the same as the input type. In this way the routine is much like many of the intrinsic functions defined as part of the Fortran standard. An example program is given below:

program main
use extension_mod
implicit none

integer::xi,yi
real::xr,yr
complex::xz,yz

xi = 2
xr = 2.0
xz = 2.0

yi = f(xi)
yr = f(xr)
yz = f(xz)

end program main

## Intrinsic functions

One can extend intrinsic functions. This is similar to overload operators.

Here we will demonstrate this by extending the sqrt function. The intrinsic function is not implemented for arguments of integer type. This is because there is no clear idea how to define the result of non integer type (e.g. ${\displaystyle {\sqrt {4}}=2}$, but how to define ${\displaystyle {\sqrt {5}}=2.236\ldots }$). We implement a method here where the result is always the nearest integer.

module sqrt_int_module

implicit none
private

! use intrinsic sqrt for data types which are not overloaded
intrinsic :: sqrt

! extend sqrt for integers
public :: sqrt
interface sqrt
procedure sqrt_int
end interface sqrt

contains
pure function sqrt_int( int ) result( res )
integer, intent(in) :: int
integer :: res

res = nint( sqrt( real(int) ) )

end function sqrt_int
end module sqrt_int_module

program main
use sqrt_int_module
implicit none

integer :: i

do i = 1, 7
print *, "i:", i, "sqrt(i):", sqrt(i)
end do

end program main

## Derived Data Types

Fortran 90 and newer supports the creation of new data types which are composites of existing types. In some ways this is similar to an array, but the components need not be all of the same type and they are referenced by name, not index. Such data types must be declared before variables of that type, and the declaration must be in scope to be used. An example of a simple 2d vector type is given below.

type::vec_t
real::x,y
end type

Variables of this type can be declared much like any other variable, including variable characteristics such are pointer or dimension.

type(vec_t)::a,b
type(vec_t),dimension(10)::vecs

Using derived data types, the Fortran language can be extended to represent more diverse types of data than those represented by the primitive types.

Operators can be overloaded so that derived data types support the standard operations, opening the possibility of extending the Fortran language to have new types which behave nearly like the native types.

### Assignment

The assignment operator = can be overloaded. We will demonstrate this by the following example. Here, we define how the assignment of a logical type on the left and an integer on the right should be performed.

module overload_assignment

implicit none
private

public :: assignment(=)
interface assignment(=)
procedure logical_gets_integer
end interface assignment(=)

contains
subroutine logical_gets_integer( log, i )
logical, intent(out) :: log
integer, intent(in) :: i

if ( i==0 ) then
log = .true.
else
log = .false.
end if
end subroutine logical_gets_integer

program main
implicit none

logical :: log

log = 0
print *, "log=0:", log          ! yields: T

log = 1
print *, "log=1:", log          ! yields: F

end program main

### Intrinsic operators

One can overload intrinsic operators, such as +,-,*.

In the following example we will overload the * operator to work as the logical .and..

module overload_operator

implicit none
private

public :: operator(*)
interface operator(*)
procedure logical_and
end interface operator(*)

contains
pure function logical_and( log1, log2 ) result( log_res)
logical, intent(in) :: log1, log2
logical :: log_res

if ( log1 .and. log2 ) then
log_res = .true.
else
log_res = .false.
end if
end function logical_and

program main
implicit none

logical, parameter :: T = .true., F = .false.

print *, "T*T:", T*T                  ! yields: T
print *, "T*F:", T*F                  ! yields: F
print *, "F*T:", F*T                  ! yields: F
print *, "F*F:", F*F                  ! yields: F

end program main

### New operators

One can create newly self-created operators.

We demonstrate this by the following example: We create an operator .even. <int> which outputs a logical if the given integer is even.

module new_operator

implicit none
private

public :: operator(.even.)
interface operator(.even.)
procedure check_even
end interface operator(.even.)

contains
pure function check_even( int ) result( log_res )
integer, intent(in) :: int
logical :: log_res

if ( modulo( int, 2 )==0 ) then
log_res = .true.
else
log_res = .false.
end if
end function check_even
end module new_operator

program main
use new_operator
implicit none

integer :: i

do i = 1, 6
print *, "i:", i, "even?", .even. i
end do

end program main


text:
Waveforms
Spring 2019

ECE 110

Course Notes

Learn It!

Required

Explore More

Learn It!

An electrical waveform is an electrical quantity (such as voltage, current or power) that can vary over time. Waveforms can deliver power. For example, in the United States, the voltage across the terminals of a power outlet is a sinusoidal waveform with frequency 60 Hz, which is equivalent to a period of 1/60 s. Waveforms can also convey information. In lab, you use a waveform called pulse width modulation (PWM) to spin a motor at a desired speed. Waveforms that transmit information and are also called signals.

Figure 1

Fig. 1: Electroencephalogram (EEG). The upper waveform is the first published electrical brain wave recording of a human (Hans Berger, 1924). The lower waveform is a 10 Hz timing sinusoid (so it has period 0.1 s). Image source.

We can analyze circuits involving time-varying waveforms because Ohm's law, KVL, KCL and the power formulas are all valid at every instant of time. For example, if a voltage waveform $V(t)$ is applied across a resistor $R$, then the resistor's current waveform $I(t)$ and power waveform $P(t)$ satisfy the following equations.
\begin{align}
I(t) &= \frac{V(t)}{R} & & \text{at each time }t \\
P(t)&=V(t)I(t)=(I(t))^2 R=\frac{(V(t))^2}{R} & & \text{at each time }t
\end{align}

We often wish to summarize the overall behavior of a waveform by a single number. For a periodic power waveform of period $T$, we calculate the average power $P_{\text{avg}}$ as follows. First integrate $P(t)$ from $0$ to $T$ to obtain the "area under the curve" in one period. Then divide the result by $T$ to give the average (or mean) value during the period.
\begin{align}
P_{\text{avg}}&=\frac{1}{T} \int_0^T P(t) dt \label{WVF-API}
\end{align}
To summarize a periodic voltage or current waveform of period $T$, we do not use the average. Instead, we use the root mean square (RMS) voltage $V_{\text{rms}}$ or current $I_{\text{rms}}$. RMS is literally the square root of the mean of the square of the waveform.
\begin{align}
V_{\text{rms}}&=\sqrt{\frac{1}{T} \int_0^T (V(t))^2 dt} \label{WVF-VRM} \\
I_{\text{rms}}&=\sqrt{\frac{1}{T} \int_0^T (I(t))^2 dt} \label{WVF-IRM}
\end{align}
The reason for this definition is that it gives the average power across a resistor $R$ similar formulas as for the nonvarying case.
\begin{align}
P_{\text{avg}}=V_{\text{rms}}I_{\text{rms}}=\frac{(V_{\text{rms}})^2}{R}=(I_{\text{rms}})^2 R \label{WVF-PFS}
\end{align}
We derive $P_{\text{avg}}=(V_{\text{rms}})^2/R$ below. The other derivations are similar.
\begin{align}
P_{\text{avg}}&=\frac{1}{T} \int_0^T P(t) dt & & \text{by the definition of } P_{\text{avg}} \text{ in } \eqref{WVF-API} \\
&= \frac{1}{T} \int_0^T \frac{(V(t))^2}{R} dt & & \text{since } P(t)=\frac{(V(t))^2}{R} \text{ at each time }t \\
&= \frac{1}{R} \frac{1}{T} \int_0^T (V(t))^2 dt \\
&= \frac{(V_{\text{rms}})^2}{R} & & \text{by the definition of } V_{\text{rms}} \text{ in } \eqref{WVF-VRM}
\end{align}
Note that, in general, $V_{\text{avg}}$ and $I_{\text{avg}}$ have no relationship with $P_{\text{avg}}$.

PWM is a technique of delivering a specific amount of average power to a motor (or other device) so that it operates at a desired rate. We show that the average power delivered by PWM to a resistor is directly proportional to the PWM's duty cycle.

Figure 2

Fig. 2: Pulse width modulation (PWM) voltage signal. This waveform $V(t)$ has amplitude $A$ volts, period $T$ seconds and duty cycle $T_{\text{on}}/T$.

We begin evaluating $V_{\text{rms}}$ by plotting $(V(t))^2$ and visually integrating it over one period.

Figure 3

Fig. 3: PWM voltage squared. This waveform $(V(t))^2$ is the square of the PWM in Fig. 2. The shaded area represents $\int_0^T (V(t))^2 dt = A^2 T_{\text{on}}$.

\begin{align}
V_{\text{rms}}&=\sqrt{\frac{1}{T} \int_0^T (V(t))^2 dt} & & \text{by the definition of } V_{\text{rms}} \text{ in } \eqref{WVF-VRM}\\
&= \sqrt{\frac{1}{T} A^2 T_{\text{on}}} & &\text{from Fig. 3} \\
&= A\sqrt{\frac{T_{\text{on}}}{T}} \label{WVF-PV2}
\end{align}
We now obtain average power $P_{\text{avg}}$ delivered to a resistor $R$ under the PWM voltage waveform $V(t)$ shown in Fig. 2.
\begin{align}
P_{\text{avg}}&= \frac{(V_{\text{rms}})^2}{R} = \left( A\sqrt{\frac{T_{\text{on}}}{T}} \right)^2 \frac{1}{R} = \frac{A^2}{R} \frac{T_{\text{on}}}{T} \label{WFV-PPW}
\end{align}
Therefore, the RMS of a PWM waveform is proportional to the square root of its duty cycle, while the average power it delivers to a resistor is directly proportional to its duty cycle.

Sinusoidal waveforms are common in electrical engineering. For example, an AC generator produces a sinusoidal voltage as a result of its rotational motion.

Figure 4

Fig. 4: Sinusoidal voltage waveform. This waveform $V(t)=A \sin(2 \pi t/T)$ has amplitude $A$ volts and period $T$ seconds.

First we apply integration to show that $V_{\text{rms}}=A/\sqrt{2}$ for a sinusoidal voltage with amplitude $A$.
\begin{align}
V_{\text{rms}}&=\sqrt{\frac{1}{T} \int_0^T A^2 \sin^2 \left( \frac{2 \pi t}{T} \right) dt} & &\text{since } V(t)=A \sin \left( \frac{2 \pi t}{T} \right) \\
&= \sqrt{\frac{1}{T} \int_0^T A^2 \frac{1}{2} \left( 1-\cos\left( 2 \frac{2 \pi t}{T} \right) \right) dt} & &\text{since } \sin^2\theta=\frac{1}{2}(1-\cos 2\theta) \\
&= \frac{A}{\sqrt{2}} \sqrt{\frac{1}{T} \int_0^T \left( 1-\cos\left( \frac{4 \pi t}{T} \right) \right) dt} \\
&= \frac{A}{\sqrt{2}} \sqrt{\frac{1}{T} \left[ t- \frac{T}{4 \pi}\sin\left( \frac{4 \pi t}{T} \right) \right]_0^T} & &\text{since } \int \cos\theta d\theta=\sin\theta+C \\
&= \frac{A}{\sqrt{2}} \sqrt{\frac{1}{T} T} & &\text{since } \sin(4\pi)=\sin(0)=0 \\
&= \frac{A}{\sqrt{2}} \label{WVF-AS2}
\end{align}
Now that we know $V_{\text{rms}}$, we can obtain $P_{\text{avg}}$ absorbed by a resistor $R$ under the applied voltage $V(t)$ shown in Fig. 4.
\begin{align}
P_{\text{avg}}&= \frac{(V_{\text{rms}})^2}{R} = \left( \frac{A}{\sqrt{2}} \right)^2 \frac{1}{R} = \frac{A^2}{2R} \label{WFV-PSN}
\end{align}
Another way to determine $P_{\text{avg}}$ directly is to derive and plot $P(t)$ and then read off the average value $P_{\text{avg}}$ graphically.
\begin{align}
P(t)&=\frac{(V(t))^2}{R} \\
&=\frac{A^2}{R} \sin^2 \left( \frac{2 \pi t}{T} \right) & &\text{since } V(t)=A \sin \left( \frac{2 \pi t}{T} \right) \\
&=\frac{A^2}{R} \frac{1}{2} \left( 1-\cos\left( 2 \frac{2 \pi t}{T} \right) \right) & &\text{since } \sin^2\theta=\frac{1}{2}(1-\cos 2\theta) \\
&=\frac{A^2}{2R} - \frac{A^2}{2R} \cos\left( \frac{4 \pi t}{T} \right) \label{WVF-EQP}
\end{align}

Figure 5

Fig. 5: Reading $P_{\text{avg}}$ from plot of $P(t)$. From equation $\eqref{WVF-EQP}$, the power $P(t)$ absorbed by $R$ consists of a constant term $A^2/(2 R)$ and a cosine term with the same value $A^2/(2 R)$ as its amplitude. So, $P_{\text{avg}}$ is equal to the constant term $A^2/(2R)$, just as in equation $\eqref{WFV-PSN}$.

The voltage across the terminals of a power outlet in the United States is referred to as 120 V. But this value is the RMS voltage of the sinusoidal waveform, so in fact the amplitude $A=120\sqrt{2}$, which is approximately 170 V.

Figure 6

Fig. 6: Outlet voltage waveform. United States outlet voltage is 120 V RMS at 60 Hz frequency. In other words, it has amplitude approximately 170 V and period 1/60 s.

Notice that $V_{\text{avg}}=0$ for this waveform and is therefore unrelated to $P_{\text{avg}}=A^2/(2R)$.

Explore More!


text:
# 3 cards are dealt from a well shuffled deck.

1.Find the chance that none of the cards are hearts.

The answer is $\frac{39}{52}$ $\cdot$ $\frac{38}{51}$ $\cdot$ $\frac{37}{50}$

However, why can't we use complement rule here:

1- p(chance that all the cards are hearts)= 1- ($\frac{13}{52}$ $\cdot$ $\frac{12}{51}$ $\cdot$ $\frac{11}{50}$)

2.Find the chance that the cards are not all hearts.

The answer given is: 1- ($\frac{13}{52}$ $\cdot$ $\frac{12}{51}$ $\cdot$ $\frac{11}{50}$)

here I'm little confused as to why the complement rule was used.

for chance that the cards are not all hearts why can't we use:

P(1 heart) + P(2 heart) + P(No hearts)= ($\frac{13}{52}$ $\cdot$ $\frac{39}{51}$ $\cdot$ $\frac{38}{50}$) + ($\frac{13}{52}$ $\cdot$ $\frac{12}{51}$ $\cdot$ $\frac{39}{50}$) + ($\frac{39}{52}$ $\cdot$ $\frac{38}{51}$ $\cdot$ $\frac{37}{50}$)

• The complement of "none of the cards are hearts" is NOT "all of the cards are hearts." – lulu Mar 22 '18 at 9:50

Find the probability that if three cards are drawn from a well-shuffled deck that none of the cards are hearts.

You asked why we could not use the complement rule here. We can. However, the complement of the set of outcomes in which none of the cards are hearts is not the set of outcomes in which all the cards are hearts. For instance, the hand $5\color{red}{\heartsuit}7\color{red}{\diamondsuit}J\clubsuit$ contains a heart, but not all of the cards in the hand are hearts. The complement of the set of outcomes in which none of the outcomes are hearts is the set of outcomes in which at least one of the cards is a heart.

Since the order of selection does not matter, this problem is best handled with combinations. The number of ways of selecting a subset of $k$ elements from a set with $n$ elements is $$\binom{n}{k} = \frac{n!}{k!(n - k)!}$$

Since $13$ of the $52$ cards in the deck are hearts, $52 - 13 = 39$ are not. Hence, a favorable hand consists of drawing $3$ of the $39$ cards that are not hearts when selecting $3$ of the $52$ cards in the deck. Hence, the desired probability is $$\Pr(\text{none of the cards is a heart}) = \frac{\dbinom{39}{3}}{\dbinom{52}{3}} = \frac{\dfrac{39!}{3!36!}}{\dfrac{52!}{3!49!}} = \frac{39!}{3!36!} \cdot \frac{3!49!}{52!} = \frac{39 \cdot 38 \cdot 37}{52 \cdot 51 \cdot 50}$$ We can also compute the probability by subtracting the probability of the complementary event from $1$. As stated above, the complementary event is the event that at least one of the selected cards is a heart. The number of ways of selecting $k$ of the $13$ hearts and $3 - k$ of the other $39$ cards in the deck is $$\binom{13}{k}\binom{39}{3 - k}$$ Hence, the probability that at least one of the three selected cards is a heart is $$\Pr(\text{at least one the cards is a heart}) = \frac{\dbinom{13}{1}\dbinom{39}{2} + \dbinom{13}{2}\dbinom{39}{1} + \dbinom{13}{3}\dbinom{39}{0}}{\dbinom{52}{3}}$$
Therefore, the probability that none of the selected cards is a heart is $$\Pr(\text{none of the cards is a heart} = 1 - \frac{\dbinom{13}{1}\dbinom{39}{2} + \dbinom{13}{2}\dbinom{39}{1} + \dbinom{13}{3}\dbinom{39}{0}}{\dbinom{52}{3}}$$

Find the probability that the cards are not all hearts.

This probability is found by subtracting the probability that all of the cards are hearts from $1$, which is $$\Pr(\text{not all of the cards are hearts}) = 1 - \frac{\dbinom{13}{3}}{\dbinom{52}{3}}$$

The complementary event is the probability that at least one of the cards is a heart, which we computed above.

You asked why this probability is not $$\frac{13}{52} \cdot \frac{39}{51} \cdot \frac{38}{50} + \frac{13}{52} \cdot \frac{12}{51} \cdot \frac{39}{50} + \frac{13}{52} \cdot \frac{12}{51} \cdot \frac{11}{50}$$ You corrected calculated the probability of selecting three hearts in three draws. However, when you calculated the probabilities of selecting exactly one heart in three draws or exactly two hearts in three draws, you did not account for the order in which the cards are drawn.

For instance, when you calculated the probability of selecting one heart in three draws, you calculated the probability of first selecting a heart and then selecting two non-hearts. However, the heart also could be drawn on the second or third draw, so the probability of selecting exactly one heart in three draws is actually $$\frac{\dbinom{13}{1}\dbinom{39}{2}}{\dbinom{52}{3}} = \frac{13}{52} \cdot \frac{39}{51} \cdot \frac{38}{50} + \frac{39}{52} \cdot \frac{13}{51} \cdot \frac{38}{50} + \frac{39}{52} \cdot \frac{38}{51} \cdot \frac{13}{50}$$ Similarly, the probability of selecting exactly two hearts in three draws is $$\frac{\dbinom{13}{2}\dbinom{39}{1}}{\dbinom{52}{3}} = \frac{13}{52} \cdot \frac{12}{51} \cdot \frac{39}{50} + \frac{13}{52} \cdot \frac{39}{51} \cdot \frac{12}{50} + \frac{39}{52} \cdot \frac{13}{51} \cdot \frac{12}{50}$$

• Thank you for the explanations, I think I'm getting it now. However I'm still little bit confused about the second part of question. you said " you did not account for the order in which the cards are drawn" why does order matter in this case? Your hand is still the same whether you get 2 cards first or last right? – jasminvvian Mar 23 '18 at 6:52
• The order in which the cards are drawn does not matter, which is why I used combinations. However, you considered sequences of outcomes. In particular, for the case of exactly one heart, there are three possible sequences: $HH^CH^C$, $H^CHH^C$, $H^CH^CH$. You calculated the probability of drawing a heart ($\frac{13}{52}$), then a non-heart ($\frac{12}{51}$), then another non-heart ($\frac{11}{50}$). Therefore, you have only accounted for one of the three sequences, which is why your answer is $1/3$ of the correct one. – N. F. Taussig Mar 23 '18 at 9:52

At first, I am solving only 1

As @lulu said, The complement of "none of the cards are hearts" is not "all of the cards are hearts, but "at least one card is a heart"
So, let us find the probability of at least one card being a heart.
Clearly it is $$13/52×39/51×38/50×^3C_1+13/52×12/51×39/50×^3C_2+13/52×12/51×11/50×^3C_3$$ Subtract this from $1$ and you will get your answer.

Now think about 2(tell me if you fail to think it yourself).


text:
# Exporting rotations of selected objects to a text file

I am using this script which I got from the internet and it works perfectly for 2.79. Please could someone help adding selected objects rotation to export to the same txt file? In degrees.... Thanks in advance.

# import the necessary modules we need
# in our case, blender's python API and python's os module
import bpy, os

# get the current selection
selection = bpy.context.selected_objects

# initialize a blank result variable
result = ""

# iterate through the selected objects
for sel in selection:
# get the current object's dimensions
dims = sel.dimensions
# write the selected object's name and dimensions to a string
result += "%s - %.03fm x %.03fm x %.03fm\n" % (sel.name, dims.x, dims.y, dims.z)

# get path to render output (usually /tmp\)
tempFolder = os.path.abspath (bpy.context.scene.render.filepath)
# make a filename
filename = os.path.join (tempFolder, "newfile.txt")
# open a file to write to
file = open(filename, "w")
# write the data to file
file.write(result)
# close the file
file.close()

• This is kind of hard to read – WhatAMesh Feb 27 at 7:06

You already outlined everything that needs to be done. Rotation can be reached at .rotation_euler:

Note that you can find this information by typing stuff to the console and hitting ctrl+space for auto complete.

The values are in radians so you need an easy way to convert them. You can just multiply the values by 57.29577951308232, or you could do from math import degrees or from math import degrees, radians so you have functions that can convert between the two for you.

Everything from math is imported in the console so you could just get the degrees of x rotation of the active object by typing degrees(bpy.context.object.rotation_euler.x) in the console to check that this indeed does work.

In a script this would look something like:

import bpy
from math import degrees

result = ''

for o in bpy.context.selected_objects:
#something else that you need to get to result
n = o.name
r = o.rotation_euler
result += "%s - X:%.3f, Y:%.3f, Z:%.3f \n" % (n, degrees(r.x), degrees(r.y), degrees(r.z))
# carry on to do whatever else
#write to file

• May be more robust using o.matrix_world.to_euler() will return the visual global rotation of the object. If the object is rotated via constraint for example o.rotation_euler will not reflect this. Remember too py console has convenience variables C = bpy.context, D = bpy.data – batFINGER Feb 27 at 13:12


text:
# Cross-line inline tcolorbox

As shown in the following example, underlining a cross-line sentence can be easily done by \ul via soul.sty. Can a cross-line sentence be tcolorboxed automatically in a similar way?

\documentclass{article}
\usepackage[showframe]{geometry}
\usepackage{soul}
\usepackage{tcolorbox}
\usepackage{parskip}
\begin{document}
In order to show the expected result, allow me to add a few words here.
\ul{This is a cross-line long sentence.} Underlining a cross-line sentence is easy.

In order to show the expected result, allow me to add a few words here.
\begin{tcolorbox}[colback=red!5!white, colframe=red!75!black, right=0pt,
rightrule=0pt, on line, hbox, sharp corners=east]
This is a cross-line
\end{tcolorbox}
\begin{tcolorbox}[colback=red!5!white, colframe=red!75!black, left=0pt,
leftrule=0pt, on line, hbox, sharp corners=west]
long sentence.
\end{tcolorbox}
Can a cross-line sentence be automatically tcolorboxed this way?
\end{document}

• Aug 17, 2021 at 17:29
• @leandriis, Thanks for the pointers. I did check related questions and answers (including these two) before posting my question. Aug 17, 2021 at 17:57
• Since you already went through the construction of the two boxes with the mirrored corners, adding \ul{..} inside is the simplest solution! Aug 17, 2021 at 18:03
• There is also tex.stackexchange.com/a/125319/89417 and tex.stackexchange.com/questions/251182/… that may come close. Both don't have the outline but that might be considered an advantage given that the lines in your example end up very far apart due to the height of the tcolorbox. Aug 17, 2021 at 20:09
• @SimonDispa, Great. I will think about your suggestion. Aug 18, 2021 at 14:11


text:
# Define coefficients for system of equations with 3 variables

1. Sep 9, 2009

### thecokeguy

1. The problem statement, all variables and given/known data
I have the following system og equations:

$\left\{\begin{array}{l} y + 2z = 1 \\ x + 4y + 3z = C \\ x + y + Bz = B \end{array}\right.$

a) Create the augmented matrix and reduce to row echelon form (I've already done this)
b) For which values of A and B is the system inconsistent?
c) For which values of A and B does the system have only one solution?
d) Determine the complete solution to the system for all values of A and B.

2. Relevant equations
The system in row echelon form:
$\begin{pmatrix} 1 & 4 & 3 & C \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & \frac{B - C + 3}{B + 3} \end{pmatrix}$

3. The attempt at a solution
Earlier when solving this kind of problems with systems with 2 variables, I solved it geometrically as lines being parallel, intersecting etc. Of course this could be solved geometrically too as planes, but as I can't imagine how to do it geometrically with >3 variables, I thought I've missed something.

If someone could point me in a direction, then I'll figure out the rest myself.

Thanks.

2. Sep 9, 2009

### njama

Ok, now go back of the system of equations. You got:

$$\left\{\begin{matrix} x+4y+3z=C\\ y+2z=1\\ z=\frac{B-C+3}{B+3} \end{matrix}\right.$$

What if B=-3 ?

3. Sep 9, 2009

### thecokeguy

As that's not possible, I'll guess the system would be inconsistent?

4. Sep 9, 2009

### njama

Yup, you're right.

Now find x,y,z solutions and again tell me for which is inconsistent and for which values its consistent.

5. Sep 9, 2009

### thecokeguy

I'm not sure what you want... Do you want me to back substitute to get expressions for y and x too, and then define values for B and C, which results in division by 0?

6. Sep 9, 2009

### njama

Yes, you need to do back substitution.

But where is A ? I can't see any A in the problem.

7. Sep 9, 2009

### thecokeguy

There isn't any A... Only B and C... Don't know why

8. Sep 9, 2009

### njama

Maybe C=A ? Anyway, do the back substitution and find the values of x,y,z in terms of B,C.

9. Sep 11, 2009

### thecokeguy

I don't see the relevance as x and y would only incorporate multiples of z, which would be the same fraction and still be inconsistent for B = -3 and thereby be consistent for B not equal to -3?

10. Sep 12, 2009

### njama

Yes, that's completely correct.

11. Sep 13, 2009

### thecokeguy

After 5 hours of searching. I finally found an example I could understand a general solution in.

$z=\frac{B-C+3}{B+3} \iff z(B+3) - B + C - 3 = 0$

With the equation above, it's easy to pick the values for B and C.

b) No-solutions: B=-3 and C!=0
c) One-solution: B!=-3
d) Complete solution: B=-3 and C=0 ... (-4+5z, 1-2z, z)

Is this completely wrong?

12. Sep 13, 2009

### njama

That's correct again.

For d) I would rather say infinite many solutions.

13. Sep 13, 2009

### thecokeguy

Just using the problems own words... As it's written by my professor and I like high grades, that's what I useually do :P

But thanks a lot... The feedback was just what I needed and much appreciated.

14. Sep 13, 2009

### njama

You're welcome!